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Improvements  in  equivalent  plate  modeling  of  aircraft  wings  are  presented.  Formula¬ 
tions  for  wing  mass,  stiffness,  and  loads  using  Classical  Plate  Theory  and  First  Order 
Shear  Deformation  Plate  Theory  are  given  in  a  general  manner  allowing  versatility  in  the 
selection  of  displacement  Ritz  polynomials.  A  new  technique  for  approximating  the  stiff¬ 
ness  of  an  array  of  spar  webs  with  the  stiffness  of  an  equivalent  sandwich  core  is  devel¬ 
oped.  A  formulation  allowing  wing  zones  modeled  with  Classical  Plate  Theory  and  First 
Order  Shear  Deformation  Plate  Theory  to  be  used  together  is  also  presented. 

Numerical  tests  were  performed  to  verify  the  validity  of  the  new  formulations.  Tests 
of  a  thick,  high  aspect  ratio  wing  showed  that  selecting  low  order  Ritz  polynomials  for  the 
linear  in-plane  displacements  of  a  symmetric  wing  can  lead  to  a  significant  reduction  in 
model  order  while  retaining  acceptable  accuracy.  Additional  tests  of  the  thick,  high  aspect 
ratio  wing  and  a  typical  supersonic  transport  wing  showed  that  an  array  of  spar  webs  may 
be  accurately  replaced  by  an  equivalent  core  leading  to  substantial  savings  in  computation 
time. 

Results  from  this  work  will  allow  complex  wing  configurations  to  be  modeled  using 
equivalent  plates  with  accuracy  comparable  to  finite  element  analysis,  yet  with  much 
greater  computational  efficiency.  The  results  add  support  for  the  use  of  equivalent  plate 
modeling  as  the  primary  wing  structural  analysis  tool  in  multi-disciplinary  design  optimi¬ 
zation  of  aircraft  structures. 
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Introduction 


Background 

In  the  early  days  of  mathematical  modeling  of  aircraft  wing  structures,  equivalent 
beam  models  were  used  to  suitably  represent  isotropic  high  aspect  ratio  wings.  With  the 
growing  pursuit  of  fast  subsonic  and  supersonic  flight  in  the  late  1940’s,  potential  wing 
configurations  took  on  a  vaiiety  of  shapes  with  low  aspect  ratio  and  thin  cross  section. 
The  new  wing  configurations,  which  became  subject  to  chordwise  bending,  could  no 
longer  be  accurately  modeled  using  classical  beam  theory.  Attempts  in  the  1950’s  to 
model  low  aspect  ratio  wings  as  wide  beams  and  equivalent  plates  were  essentially  aban¬ 
doned  with  the  development  of  high  speed  computing  machines  and  matrix  methods  of 
sti-uctural  analysis  (GR61).  Finite  element  (FE)  techniques  viably  emerged  from  the  new 
computing  capability.  It  made  possible  the  analysis  of  complex  low  aspect  ratio  wing 
structures  giving  stress,  deflection,  and  vibration  infoiTnation  with  acceptable  accuracy. 
The  finite  element  method  grew  in  acceptance  as  its  reliability  was  tested,  and  it  has  now 
become  industry’s  standard  analysis  tool  for  wing  design  (SCH81). 

Difficulties  were  encountered,  however,  in  the  use  of  the  finite  element  method  during 
the  preliminary  stages  of  wing  design  when  structural  shape  variations  are  examined. 
Finite  element  models  lead  to  large  matrix  equations  that  were  computationally  intensive 
even  with  the  new  computing  power.  Since  the  finite  element  mesh  had  to  change  corre¬ 
sponding  to  shape  changes,  comparison  of  many  alternative  planforms  proved  to  be  cum¬ 
bersome  and  time  consuming.  More  complexity  was  added  to  the  modeling  problem 
when  fiber  composite  materials  came  into  use.  Thus,  interest  in  the  development  of  alter¬ 
native,  computationally  efficient  wing  models,  tailored  towards  wing  shape  optimization, 
was  rekindled. 

Two  different  approaches  to  the  problem  have  emerged;  use  of  expedited  finite  ele¬ 
ment  modeling  and  use  of  equivalent  plate  modeling.  Major  advances  have  been  made  in 
improving  finite  element  modeling  for  use  in  design  oriented  structural  analysis  (HG92). 
In  the  area  of  airframe  structures,  numerical  and  analytical  methods  were  used  to  obtain 
design  sensitivity  information  with  respect  to  both  sizing  and  shape  design  variables 
(HA93,LAN94).  Another  improvement  has  been  the  development  of  integrated  FE  /  auto- 
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mated  mesh  generation,  which  allows  changes  in  planform  and  airfoil  shapes  to  be  han¬ 
dled  efficiently  (HA93).  Although  advances  in  modeling  and  improvements  in 
computational  power  have  made  practical  the  shape  optimization  of  finite  element  models 
having  many  design  variables,  the  computational  price  remains  high.  In  multi-disciplin¬ 
ary  design  optimization,  where  the  structural  analysis  is  just  part  of  a  wider  optimization 
scheme,  the  cost  of  using  a  finite  element  approach  may  still  be  too  expensive  (LIV90). 

The  full  capability  of  wing  structural  modeling  based  on  equivalent  plates  is  still 
being  explored.  Equivalent  plate  models  have  been  found  capable  of  capturing  spanwise 
and  chordwise  bending  of  various  planforms  using  a  small  number  of  degrees  of  freedom. 
An  equivalent  plate  model  based  on  Classical  Plate  Theory  (CPT)  was  used  in  the 
aeroelastic  tailoring  and  structural  optimization  (TSO)  computer  program 
(LRB77,MC83).  This  simplified  capability  has  been  used  extensively  for  preliminary 
wing  design  with  good  results  compared  to  finite  element  analysis  for  real  wings  such  as 
those  of  the  YF16  and  F15  fighter  aircraft  (LRB77,TR80).  The  high  computational  effi¬ 
ciency  of  TSO  coupled  with  its  limitation  of  handling  only  a  single  trapezoidal  planform 
motivated  the  development  of  an  equivalent  plate  approach  capable  of  handling  planforms 
composed  of  multiple  trapezoidal  sections,  with  and  without  symmetry  about  the  wing 
mid-surface  (GI86,GI89).  Better  modeling  of  the  internal  wing  structure  was  also 
included.  Results  showed  that  equivalent  plate  models  could  accurately  predict  stresses  as 
well  as  deflections  and  mode  shapes  for  wings  with  complex  planforms  while  offering  sig¬ 
nificant  savings  in  computer  time.  This  capability  was  extended  to  include  analytically 
derived  design  sensitivities  with  respect  to  the  shape  and  sizing  design  variables  (LIV90). 
It  was  also  shown  that  the  planform  shape  could  be  modeled  such  that  all  integration  could 
be  done  analytically,  thus  eliminating  the  need  for  coordinate  transformations  and  numeri¬ 
cal  integration. 

A  study  of  equivalent  plate  modeling  applied  to  a  potential  High  Speed  Civil  Trans¬ 
port  (HSCT)  wing  revealed  some  limitations  of  plate  models  based  on  CPT  (LSB93). 
Results  from  the  study  showed  that  a  CPT  based  equivalent  plate  model  does  not  give 
good  results  for  a  wing  with  few  spars  and  ribs  and  low  transverse  shear  stiffness.  In  light 
of  this  finding,  a  new  equivalent  plate  modeling  capability  was  developed,  based  on  First 
Order  Shear  Deformation  Plate  Theory  (FSDPT),  which  showed  much  better  correlation 
with  finite  element  results  for  the  HSCT  wing  (LIV194,LIV294,LSB93). 
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Motivation 

The  newly  developed  equivalent  plate  structural  modeling  approach  using  FSDPT 
provided  more  accurate  results  than  a  CPT  based  model  for  certain  wing  configurations. 
Yet  in  using  FSDPT,  some  of  the  computational  efficiency  sought  in  using  a  plate  model¬ 
ing  approach  was  lost.  The  FSDPT  capability  (in  its  form  presented  in  LTV  194)  had  limi¬ 
tations  that  led  to  increased  computational  cost.  Among  the  causes  for  this  increased 
computational  effort  was  the  way  in  which  polynomial  series  were  chosen  to  represent 
displacement  functions  in  the  model.  Originally,  all  displacement  functions  had  to  be  rep¬ 
resented  by  x-y  polynomials  of  the  same  order  (LIV194).  Since  five  displacement  fields 
are  involved  with  FSDPT  modeling  compared  with  one  field  (w  deflection)  in  CPT,  this 
resulted  in  a  significant  increase  of  the  size  of  the  system  of  equations.  Also,  a  significant 
computational  burden  was  carried  by  the  original  FSDPT  approach  because  of  the  inclu¬ 
sion  of  spar  webs  and  rib  webs  in  the  model.  The  computational  time  required  to  calculate 
the  mass  and  stiffness  contributions  from  an  array  of  many  webs  became  very  time  con¬ 
suming. 

For  a  low  aspect  ratio  wing,  like  that  of  the  HSCT,  which  has  a  complex  planform 
geometry,  low  transverse  sheai'  stiffness,  and  multiple  control  surfaces,  equivalent  plate 
models  based  on  CPT  alone  and  FSDPT  alone  ai'e  inadequate.  It  has  already  been  shown 
that  a  CPT  based  model  cannot  capture  the  transverse  shear  affects  of  such  a  wing 
(LSB93).  A  FSDPT  based  model  gave  good  results  with  added  computational  cost.  How¬ 
ever,  the  computational  cost  would  grow  significantly  if  several  control  surfaces  also  mod¬ 
eled  using  FSDPT  were  included. 

Facing  these  limitations  of  the  current  equivalent  plate  capabilities,  it  is  important  to 
make  improvements  in  plate  modeling  to  retain  the  accuracy  of  the  FSDPT  approach 
while  bolstering  computational  efficiency.  Specifically,  it  is  sought  to  generalize  the  cur¬ 
rent  FSDPT  capability  to  handle  multiple  “zones”  (LSB93)  and  allow  the  order  of  each 
displacement  polynomial  to  be  independently  assigned.  To  further  reduce  computation 
time  it  is  desired  to  approximate  the  stiffness  contribution  of  an  array  of  spar  and  rib  webs 
with  the  stiffness  of  an  equivalent  sandwich  core  (LIV194).  Thus,  instead  of  calculating 
stiffness  and  mass  contributions  of  spar  and  rib  webs  one  by  one,  a  single  equivalent  core 
is  evaluated  once  per  new  wing  shape.  Finally  sought  is  a  new  capability  allowing  zones 
based  on  FSDPT  to  be  used  in  conjunction  with  zones  based  on  CPT,  thus  allowing  wings 
such  as  that  of  the  HSCT  with  control  surfaces  to  be  accurately  and  efficiently  modeled 
using  equivalent  plates. 
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Scope  of  Work 

Following  this  introduction,  Chapter  1  discusses  the  conventions  and  Ritz  solution 
method  used  in  the  wing  modeling  problem.  Chapters  2  and  3  present  the  equivalent  plate 
formulations  for  wing  structural  mass,  structural  stiffness,  and  loading  conditions  using 
Classical  Plate  Theory  and  First  Order  Shear  Deformation  Plate  Theory  respectively. 
Both  formulations  are  derived  for  numerical  analysis  such  that  they  allow  the  order  of 
each  displacement  polynomial  to  be  independently  chosen.  Chapter  4  discusses  zone  con¬ 
nections  and  boundary  conditions,  including  a  new  capability  for  combining  CPT  and 
FSDPT  zones.  Chapter  5  presents  a  method  used  to  approximate  the  stiffness  of  a  spar 
web  array  by  the  stiffness  of  a  sandwich  core.  Chapter  6  compares  the  results  from  use  of 
the  new  plate  modeling  capabilities  with  available  data  for  two  different  test  wings.  Chap¬ 
ter  7  concludes  with  a  summary  and  discussion  of  possible  extensions  of  this  work. 

Detailed  derivations  of  the  equivalent  plate  formulations  in  Chapters  2  and  3  may  be 
found  in  Appendices  A  and  B.  Appendix  C  contains  analytical  expressions  for  all  inte¬ 
grals  encountered  in  the  wing  math  model.  Appendix  D  contains  data  and  figures  from  the 
test  wings.  Appendix  E  contains  information  concerning  the  use  of  the  program  CON¬ 
NECT  and  supporting  subroutines  developed  during  this  research  effort. 


Chapter  1 : 
Wing  Model 


1.1  Overview 

In  this  chapter  the  paiameters  for  a  mathematical  wing  model  based  on  equivalent 
plate  theory  are  established,  and  the  solution  technique  is  discussed.  The  chapter  begins 
by  discussing  the  configuration  of  a  general  wing  box  structure  and  defining  the  conven¬ 
tion  used  for  the  model’s  shape  design  variables  based  on  planform  geometry  and  wing 
depth  distribution.  Next  discussed  is  the  convention  used  to  express  wing  constiuction 
parameters  as  polynomial  series,  where  the  coefficients  are  the  model  s  sizing  design  vari¬ 
ables.  Wing  material  properties  are  discussed,  and  constitutive  relationships  are  given  for 
the  wing  structural  components.  Following  this  is  a  presentation  of  the  Ritz  solution  tech¬ 
nique  as  applied  to  the  wing  model.  Here  the  stiffness  matrix  niass  matrix  [AT],  gener¬ 
alized  displacements  vector  {g},  and  generalized  load  vector  {P}  are  defined.  Briefly 
discussed  are  details  of  the  static  solution  and  eigensolution.  The  chapter  concludes  with 
a  discussion  of  output  grid  specification. 

1.2  General  Wing  Box  Configuration 

The  general  wing  box  configuration  considered  in  this  work  consists  of  skins,  spar 
caps,  rib  caps,  spar  webs,  and  rib  webs  as  shown  in  Figure  1.1.  Skins  may  be  constructed 
of  multiple  unidirectional  fiber  composite  laminae.  The  thickness  and  distance  from  the 
wing  mid-plane  of  each  skin  lamina  may  vary  over  the  area  of  wing.  Stiffening  spars  and 
ribs  are  attached  to  the  upper  and  lower  skins  generally  in  some  systematic  array.  Along 
the  length  of  a  spar  or  rib  the  cap  distance  from  the  wing  mid-plane  may  vary.  The  cap 
areas  of  both  spars  and  ribs  also  may  vary  linearly  along  their  lengths.  Spar  webs  and  rib 
webs,  like  skins,  may  also  be  constructed  of  multiple  composite  layers  and  are  intended 
primaiily  to  provide  shear  stiffness.  Spar  webs  connect  spar  caps  on  the  upper  surface 
with  their  parallels  on  the  lower  surface.  Rib  caps  do  the  same  for  parallel  rib  caps. 
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Figure  1.1:  Wing  Box  Configuration 

1 .3  Wing  Planform  Geometry 

The  geometry  of  a  general  wing  planform  is  defined  in  a  global  x-y-z  coordinate  sys¬ 
tem.  The  x-y  plane  (z=0)  corresponds  to  the  reference  surface  of  the  wing.  A  typical  wing 
planform  is  shown  in  Figure  1.2.  The  whole  planform  may  be  divided  into  several 
“zones”.  Each  zone  may  have  its  geometry  defined  in  its  own  zonal  coordinate  system 
and  may  use  unique  Ritz  polynomial  functions  to  represent  its  structural  displacements. 
All  zones,  however,  share  a  common  z  axis. 


zone  2 


Figure  1.2:  Multi-zone  Wing  Pianform 
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Figure  1.3:  Skin  Panei  Shape  Design  Variabies 

Each  zone  is  composed  of  one  or  more  upper  and  lower  trapezoidal  panels.  Each 
panel  is  defined  by  6  shape  design  variables  yi,  yn,  Xpi,  XpR,  x^i,  and  x^r  as  shown  in  Figure 
1 .3.  The  vaiiables  yi  and  y^  define  respectively  the  left  and  right  y  coordinates  of  the  sides 
of  the  panel  which  are  parallel  to  the  x-axis.  The  variables  Xpi  and  Xpp  define  respectively 
the  front  left  and  front  right  x  coordinates  of  the  panel.  Likewise,  the  variables  x^p  and  x^r 
define  respectively  the  aft  left  and  aft  right  jc  coordinates  of  the  panel.  A  zone  generally 
will  also  have  associated  spars,  ribs,  and  webs.  Each  spar  is  defined  by  the  endpoints  of  a 
line  segment.  The  4  shape  design  variables  sjp,  sy^,  sxp,  and  sxr  define  a  spar  line’s  geom¬ 
etry  as  shown  in  Figure  1.4.  A  spar  is  generally  at  some  angle  A  to  the  y-axis.  Ribs  are 


y 


(rxAoryi) 


Figure  1.4:  Spar  and  Rib  Shape  Design  Variables 
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parallel  to  the  x-axis  and  are  defined  by  7  shape  design  variables  ry^^,  ry^,  rxpi,  rxpR,  ry^p, 
ry^p,  and  ypjp.  The  first  6  of  these  variables  define  a  trapezoid  in  the  same  manner  as  the 
panel  ttapezoids  are  defined.  The  variable  y^^gthen  uniquely  defines  a  rib  line  of  constant 
y  value,  whose  ends  lie  on  the  front  and  aft  lines  of  the  trapezoid  as  shown  in  Figure  1.4. 
Spar  and  rib  web  lines  match  the  lines  of  the  spar  and  rib  caps  they  connect. 

Zones  may  also  have  associated  “attach  lines”,  “attach  points”,  and  “free  moving 
points”.  An  attach  line  identifies  the  common  surface  between  two  adjoining  zones.  In  a 
particular  zone’s  coordinate  system,  an  attach  line  is  defined  by  the  4  variables  atyp,  aty^, 
atXp,  and  atXp  which  give  the  left  and  right  x  and  y  coordinates  of  the  line’s  endpoints.  An 
attach  point  lying  on  an  attach  line  is  defined  simply  by  its  x-y  coordinates.  At  attach 
points  concentrated  computational  springs  may  be  used  to  impose  boundary  conditions  or 
to  represent  the  local  stiffness  of  an  attachment  between  zone  boundaries  (LIV90, 
LIV194).  Free  moving  points  are  also  defined  by  their  x-y  coordinates  and  are  used  to  rep¬ 
resent  points  on  the  wing  where  mass  is  lumped  or  where  concentrated  force  loads  are 
applied  (LIV90). 


1.4  Wing  Depth  and  Construction 


The  depth  distribution  of  a  wing  component  is  defined  as  its  distance  from  the  z=0 
reference  plane.  A  depth  distribution,  expressed  as  a  simple  polynomial  series,  may  be 
defined  across  an  entire  zone  for  each  skin  layer  and  each  array  of  spars  and  ribs.  For 
example,  the  depth  distribution  across  a  single  zone  for  the  jlth  skin  layer  is  given  by 


Nhj, 


hji(x,y)  = 


mhi^  nhj^ 


k  -  1 


(1.1) 


where  the  coefficients  H-,  are  shape  design  variables,  and  the  powers  mh^.  and  nh,,  may  be 
chosen  to  specify  a  polynomial  in  x  and  y  containing  Nhji  terms  (Actually,  the  powers  asso¬ 
ciated  with  the  ^h  term  of  the  yVth  layer  should  be  denoted  rnkj^  ,  nkj^  ).  For  a  wing  whose 
skin  thickness  is  small  compared  to  its  depth  it  is  justifiable  to  use  only  2  depth  distribu¬ 
tions  per  zone.  In  such  a  case,  all  skins,  spars,  and  ribs  on  the  upper  surface  may  be 
assigned  one  distiibution  hu(x,y)  while  those  on  the  lower  surface  may  have  a  second  dis¬ 
tribution  hi^ix,y)  (LrV194). 

Upper  and  lower  wing  skins  may  be  composed  of  several  layers  of  fiber  composite 
material.  Each  layer  may  have  its  principal  fiber  direction  oriented  at  some  angle  p  to  the 
X  axis.  For  each  trapezoidal  panel  the  thickness  of  the  jlth  layer  may  be  defined  by  a  sim- 
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pie  polynomial  series  composed  of  Ntji  terms  and  having  the  form 


where  the  coefficients  fji^  are  sizing  design  variables,  and  the  powers  w4  and  «4  are  pre¬ 
assigned.  Again  it  should  be  mentioned  that  the  powers  associated  with  term  k  in  layer  yV 
should  be  denoted  mt-,  ,nt;  .  However,  the  simpler  notation  in  Eq.  1.2  is  adequate  at  this 

J^k  J^k 

point  for  illustrative  purposes. 

Spar  and  rib  cap  areas  may  vary  linearly  along  their  respective  longitudinal  axes.  The 
cap  area  of  the  jsth  spar  may  be  expressed  as  a  linear  function  of  y  by 


(y) 


(1.3) 


where  the  coefficients  and  A^^j  are  sizing  design  variables.  The  cap  area  of  the 

7>th  rib  varies  linearly  with  x  and  may  be  expressed  as 


(x) 


(1.4) 


where  the  coefficients  A^^q  and  A^^^  are  sizing  design  variables. 

Spar  and  rib  webs  are  constructed  similarly  to  the  wings  skins  in  that  each  composite 
layer  of  a  web  may  have  a  varying  thickness  and  independently  oriented  fiber  direction. 
Equation  1.2  may  be  used  to  express  the  thickness  of  a  single  web  layer.  However,  for 
purposes  of  this  work  it  is  desirable  to  treat  the  thickness  of  each  layer  as  varying  linearly 
like  the  cap  areas  and  having  just  2  polynomial  coefficients.  The  spar  web  thickness  then 
varies  linearly  in  y  while  the  rib  web  thickness  varies  linearly  in  x. 

Other  sizing  design  variables  include  concentrated  masses  and  the  stiffness  coeffi¬ 
cients  of  concentrated  springs.  A  concentrated  mass  represents  a  lumped  mass  applied  at 
a  free  moving  point  or  an  attach  point  as  discussed  in  the  previous  section.  Concentrated 
springs  may  be  used  to  represented  local  stiffness  at  an  attach  point.  Large  spring  stiffness 
coefficients  may  be  used  in  a  penalty  method  to  impose  a  zero  displacement  and  zero  rota¬ 
tion  boundary  condition  or  to  enforce  continuity  at  a  point  common  to  2  zones.  Rotational 
springs  may  be  used  to  represent  hinge  and  actuator  stiffness.  Carefully  selected  spring 
stiffness  coefficients  may  be  used  to  compensate  for  differences  between  numerical  testing 
and  experiment  by  representing  the  actual  flexibility  of  an  experimental  support  setup 
(LIV 1 94,LRB77,LSF88,KA92,GR6 1). 
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1 .5  Material  Properties 

The  material  density  p  of  each  wing  component  is  used  in  determining  the  wing’s  the¬ 
oretical  mass.  Material  stiffness  properties  are  used  in  constitutive  laws  for  each  wing 
component  to  determine  the  wing’s  stiffness. 

The  constitutive  relationship  between  in-plane  stresses  <J  and  strains  e  for  they7th  skin 
layer  is  expressed  as 


a 

8 

XX 

XX 

a 

yy 

^  =  [2]/ 

8 

yy 

.  a 

xy 

^  yxy  ^ 

where  the  material  matrix  [Qh  contains  the  stiffness  properties  of  the  composite  layer 
appropriately  transformed  into  the  x-y  axes  system  based  on  the  orientation  of  the  layer’s 
principle  material  axes.  Each  spar  cap  and  rib  cap  is  treated  as  having  uniaxial  stiffness 
based  the  cap’s  equivalent  Young’s  modulus  E. 

A  spar  web  lies  in  the  ri-z  plane  where  Tj  is  the  coordinate  of  the  web’s  associated 
spar  line  as  seen  in  Figure  1.5.  As  discussed  in  Chapter  3,  it  is  assumed  that  e^=0.  There¬ 
fore,  the  constitutive  relationship  for  the  jlth  web  layer  is  expressed  as 


(1.6) 


where  the  material  matrix  [U\ji  contains  the  stiffness  properties  of  the  composite  layer 


Figure  1.5:  Spar  Web  and  its  Associate  Axes 
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appropriately  transformed  into  the  Ti-z  coordinate  system  based  on  the  orientation  of  the 
layer’s  principle  fiber  dkection  (LrV194).  A  rib  web  is  treated  in  a  similar  manner  except 
that  it  lies  in  an  x-z  plane.  The  constitutive  relationship  for  the  jlth  layer  of  a  rib  web  is 
given  by 


(1.7) 


The  stress  in  a  web  can  be  retrieved  once  the  deformation  and  strains  are  determined. 


1.6  Ritz  Solution  Technique 

The  Ritz  method  is  used  to  obtain  an  approximate  solution  for  the  generalized  dis¬ 
placements  that  minimize  the  total  energy  of  the  wing  box  structure.  The  total  energy 
associated  with  the  wing  model  is 


where  U  is  the  potential  energy  stored  in  the  structure  through  deformation,  Q  is  the  work 
of  the  applied  loads  moving  through  the  corresponding  structural  deflections,  and  T  is  the 
kinetic  energy  associated  with  the  mass  of  the  structure  (GI86).  These  energy  terms  are 
expressed  as  functions  of  the  displacements  of  the  structure.  All  structural  displacements 
can  be  individually  expressed  as  a  polynomial  series  of  the  general  form 


Nc 

_  me-  nci 

f{x,y,t)  =  r ^  y 

i  =  1 


(1.9) 


where  the  coefficients  c{t\  are  unknown  and  vary  with  time,  and  the  powers  mt,  and  nt,  are 
predetermined  from  possible  polynomial  combinations  of  increasing  order  as  shown  in 
Figure  1.5.  Since  the  displacement  polynomials  are  used  in  determining  the  energy  of  the 
structure,  the  energy  terms  U,  Q,  T,  and  thus  become  functions  of  the  unknown  c(0r 
The  Ritz  solution  requires 


(1.10) 
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Figure  1.6:  Complete  Polynomials  Through  7th  Order 


Application  of  Eq.  1.10  to  a  vibrating  wing  under  harmonic  excitation  results  in  a  system 
of  simultaneous  equations  expressed  in  matrix  form  as 


where  is  the  symmetric  stiffness  matrix  derived  from  the  potential  energy,  [M\  is  the 
symmetric  mass  matrix  derived  from  the  kinetic  energy,  co  is  the  structure’s  natural  vibra¬ 
tional  frequency,  {P}  is  the  generalized  load  vector  derived  from  the  applied  load  condi¬ 
tions,  and  {c}is  the  vector  of  unknown  generalized  displacements  (GI86).  Chapters  2  and 
3  discuss  how  the  stiffness  matrix,  mass  matrix,  and  load  vector  are  derived  for  a  particu¬ 
lar  wing  using  different  plate  theories. 


1.7  Static  Solution  for  Displacements  and  Stresses 

For  a  single  static  loading  condition,  Eq.  1.11  becomes 

[k]{9}  =  {P} 


(1.12) 


where  {<7}  is  now  the  vector  of  unknown  generalized  displacements.  This  equation  holds 
for  linear  structural  systems.  The  solution  subroutines  employed  in  this  work  make  use  of 


13 


the  symmetry  and  inherent  sparsity  of  the  stiffness  matrix  [A^  (FE75).  The  solution 
method  is  based  on  a  modified  Cholesky  decomposition  technique  which  factorizes  [K\  as 


[AT]  =  [L]  [D]  [if 

(1.13) 

where  [L]  is  a  lower  triangular  matrix,  and  [D]  is  a  nonsingular  diagonal  matrix.  The  solu¬ 
tion  vector  {q}  is  then  obtained  through  a  series  of  forward  and  backward  substitutions 
using  the  factored  matrices  (FE75,LIV90). 


1.8  Eigensolution  for  Natural  Vibration  Frequencies  and  Modes 

For  a  freely  oscillating  undamped  linear  structure,  Eq.  1.11  becomes 

[a:-(d^m]{?}  =  [0] 


(1.14) 


where  {q}  is  now  the  vector  of  unknown  generalized  displacements  defining  the  vibration 
mode  shape  which  corresponds  to  the  natural  frequency  co.  The  QZ  algorithm  (GVL89)  is 
used  here  to  find  eigenvalues  and  eigenvectors  of  the  linear  generalized  eigenvalue  prob¬ 
lem  in  Eq.  1.14. 


Figure  1.7:  Pane!  Output  Grid 
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1.9  Output 

From  the  static  and  dynamic  solutions,  it  is  possible  to  give  the  static  deflection,  static 
stress,  and  modal  deflection  at  any  point  on  the  wing  model.  It  is  advantageous  to  be  able 
to  establish  a  grid  of  output  points  on  panels,  spars,  and  ribs  as  desired.  As  shown  in  Fig¬ 
ure  1.6  a  panel  output  grid  may  be  specified  over  a  square  where  -1<  ^  <  1  and  -1<  Tj  <  1 
by  choosing  the  number  of  points  desired  in  the  x  and  y  directions.  This  grid  of  points  is 
then  mapped  onto  a  selected  trapezoidal  output  panel  yielding  a  grid  of  output  points  in 
the  jc-y  axis  system.  The  same  thing  may  be  done  in  one  dimension  for  spars  and  ribs. 


Chapter  2: 

Equivalent  Plate  Formulation  Using 
Classical  Plate  Theory 


2.1  Overview 

This  chapter  contains  the  Classical  Plate  Theory  based  formulation  of  the  mass,  stiff¬ 
ness,  and  load  contributions  needed  to  solve  for  the  generalized  displacements  in  Eq.  1.1 1. 
The  chapter  begins  with  a  discussion  of  the  assumptions  of  Classical  Plate  Theory  and  the 
ensuing  displacement  and  strain  relationships.  Following  this  foundation  is  the  derivation 
of  the  wing  skin  contribution  to  the  mass  and  stiffness  of  the  system.  Next  discussed  is  the 
derivation  of  the  spar  cap  contribution  to  mass  and  stiffness.  The  derivation  of  the  rib  cap 
contribution  to  mass  and  stiffness  is  given  as  well.  Also  included  is  a  discussion  of  the 
contribution  of  concentrated  masses  to  the  mass  of  the  system.  Following  this  is  the  for¬ 
mulation  of  the  load  vector  including  concentrated  force  load  and  distributed  pressure  load 
conditions.  Finally  presented  is  the  method  used  to  calculate  displacement  and  stress  out¬ 
put  from  the  generalized  displacement  solution.  Appendix  A  includes  detailed  derivation 
of  mass  and  stiffness  contributions  from  the  wing  structural  components. 


2.2  CPT  Displacements  and  Strains 

Classical  Plate  Theory  (CPT)  is  the  name  given  to  the  small-deflection  theory  of 
bending  of  elastic  thin  plates  (RE84).  The  central  assumption  of  CPT,  known  as  the  Kir- 
choff-Love  assumption,  is  that  planes  normal  to  a  plate’s  mid-surface  will  remain  plane 
and  normal  to  the  plate’s  mid-surface  after  the  plate  is  deformed  in  bending  (RE84,J075). 
This  assumption  is  equivalent  to  ignoring  the  shear  strains  and  for  a  plate  in  the  x-y 
plane.  Because  the  plate  is  thin,  the  transverse  normal  strain  is  also  negligible  making 
this  in  essence  a  plane  stress  problem.  Because  wing  symmetry  about  the  x-y  plane  is 
assumed  here,  there  are  no  in-plane  displacements  at  points  on  the  plate’s  mid-surface. 

According  to  the  above  assumptions,  the  transverse  displacement  w  is  a  function  of  x 
and  y  only.  The  Kirchoff-Love  assumption  allows  the  displacement  u  in  the  x  direction 
and  the  displacement  v  in  the  y  direction  to  be  given  in  terms  of  w  by 
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dw 

9w 

V  =  =  -zw, 

dy  y 


(2.1) 

(2.2) 


The  in-plane  strains  e„,  Eyy,  and  are  given  by 


du 

^xx  -  -  ^^’x 


(2.3) 


dv 

Sj.  =  55;  =  -^“'■yy 


(2.4) 


du  dv  o 
~  ^  dx  ~  ~ 


(2.5) 


(LIV90,RE84).  It  is  evident  that  all  the  displacements  and  strains  are  functions  of  the 
independent  transverse  displacement  w. 

In  the  polynomial  Ritz  solution  method  described  in  Chapter  1,  it  is  necessary  to 
define  the  displacement  function  w(x,y,t)  as  a  polynomial  series  of  the  form  shown  in  Eq. 
1.9  and  having  Nq^  terms.  This  polynomial  function  may  be  expressed  in  vector  form  by 

T 

w(x,y,t)  =  {a  (x,y)}  {q^(t)} 

(2.6) 

where 

(2.7) 

and  {q^^it)}  is  a  column  vector  containing  the  time  dependent  coefficients  that  ultimately 
become  the  unknown  generalized  displacements.  Both  vectors  contain  Nq^  terms.  This 
vector  form  for  w{x,y,t)  will  be  used  in  the  mass  and  stiffness  matrix  formulations. 
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2.3  Skin  Contribution  to  Mass  and  Stiffness  for  a  CPT  Zone 

The  contribution  to  the  kinetic  energy  T  of  an  infinitesimal  skin  element  dx  by  dy  of 
the  jlth  skin  layer  of  a  single  trapezoidal  panel  is 

1  2 

dTji  =  (x,y,t)dxdy 


where  pji  is  the  constant  material  density  of  the  skin  layer,  and  tji(x,y)  is  the  thickness  of  the 
skin  layer.  Using  Eq.  2.6  for  w{x,y,t)  and  integrating  over  the  whole  trapezoidal  panel 
area  gives 

T  T 

{4jdxdy 

Summing  the  kinetic  energy  contributions  of  the  Nlu  layers  of  the  panel’s  upper  skin,  mul¬ 
tiplying  by  2  because  of  wing  symmetry,  and  summing  over  the  Np  trapezoidal  panels  in  a 
zone  leads  to  the  total  skin  mass  matrix,  x  in  dimension,  given  by 


tot 


Np  Nly 

= 

ip=ljl=l  yx 


(2.10) 


where  tji(x,y)  is  a  polynomial  function  of  x  and  y  as  given  in  Eq.  1.2.  Each  element  of 
i^skltoi  is,  thus,  a  linear  combination  of  trapezoidal  area  integrals  of  the  form 


Ij^{m,n)  =^^x^yi^dxdy 

yx 


(2.11) 


Appendix  C  discusses  how  an  integral  table  for  different  m  and  n  combinations  may  be 
analytically  constructed  for  each  trapezoidal  panel. 

The  contribution  to  the  potential  or  strain  energy  [/  of  an  infinitesimal  skin  element  dx 
by  dy  of  the  j7th  skin  layer  of  a  single  trapezoidal  panel  is 


dU-i  =  \{KV{D].i{yi}dxdy 


(2.12) 


where 


{K}  =  - 


W, 


XX 


yy 

1  2w, 


(2.13) 
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using  Eqs.  2.3,  2.4,  and  2.5;  and  [D]j,  is  the  plate  bending  stiffness  matrix  for  layer  jl 
defined  by 

(J075,  LrV90).  The  matrix  [Q]ji,  3  x  3  in  dimension,  is  the  jlth  layer’s  constitutive  matrix 
referenced  to  the  x-y  axes  as  discussed  in  Section  1.5.  Since  CPT  zones  assume  wing 
symmetry,  the  depth  distribution  h  will  be  treated  here  as  the  distance  between  the  upper 
and  lower  skins.  Referring  to  Figure  A.2,  the  jlth  skin  layer  lies  between  the  coordinates 
z=h/2-tji/2  and  z=hl2+tji/2  such  that  the  z  integral  from  Eq.  2.14  is  given  by 


rih  +  tj,)/2 


2  h 

^  =  hi 


(2.15) 


Since  it  is  assumed  that  the  wing  thickness  is  much  smaller  than  the  wing  depth  {tj/h  «1), 
the  right  hand  side  of  Eq.  2. 15  can  be  simplified  to  (LIV90).  Substituting  Eqs.  2. 13 
and  2.15  into  Eq.  2.12  and  integrating  over  the  panel  area  gives 


{gyiW]^lQ]  [W]  {qjdxdy 

vx  jl 


(2.16) 


where 


[W] 


T 


{^w’yy^ 

T 


mq.-2  nq. 

l)x  y  , 

mq.  nq.-2 

...,ng.(n^.- l)x  y 

mq^  -  1  nq^  -  1 

...,2mq.nq^x  y 


(2.17) 


(LIV90).  [W]  is  3  X  Nq^  in  dimension.  Summing  the  potential  energy  contributions  of  the 
Nly  layers  of  the  panel’s  upper  skin,  multiplying  by  2  because  of  wing  symmetry,  and 
summing  over  the  Np  panels  in  the  zone  leads  to  the  total  skin  stiffness  matrix,  Nq„  x  Nq„ 
in  dimension,  given  by 
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Np  Nly 

=  js  I  IK  fey)';/ <^•5’)  lW]^lQhdW]dxdy 

ip=ljl=lyx  K^-io) 

where  h(x,y)  and  tji(x,y)  are  polynomial  functions  of  x  and  y  as  defined  in  Eqs.  1.1  and  1.2. 
Each  element  of  [K^klut  is  a  linear  combination  of  area  integrals  of  the  form  Ij^{m,n)  as 
defined  in  Eq.  2.11. 


2.4  Spar  Cap  Contribution  to  Mass  and  Stiffness  for  a  CPT  Zone 

The  contribution  to  the  kinetic  energy  T  of  an  element  of  length  dr\  of  theyisth  spar  is 

P  2  (2.19) 

where  T)  is  the  coordinate  along  the  spar  axis  rotated  from  the  y  axis  by  an  angle  A,  is 

the  constant  material  density  of  the  spar,  and  the  cap  area  A  is  expressed  as  a  function 

xjs 

of  ri  (LIV90).  Referring  to  Figure  A.3,  all  rj  dependence  in  this  equation  may  be  changed 
to  y  dependence  using 

dr[  =  dy  - - - -  =  — ^ 

syR-syp  cos  A  (2.20) 


The  spar  cap  area  may  be  expressed  as  a  linear  function  of  y,  and  the  variable  x  may  also 
be  expressed  in  terms  of  y  using  the  spar  line  equation 


X  (y)  =  51y  +  52 


(2.21) 


Substituting  Eqs.  2.6  and  2.20  into  Eq.2.19  and  integrating  over  the  length  of  the  spar 
gives 


2cosAJy=*)'z. 


{qj 


T 

{aj  {aj 


\qjdy 


(2.22) 


Summing  the  kinetic  energy  contributions  of  the  Nsy  spars  on  the  upper  wing  surface  and 
multiplying  by  2  because  of  wing  symmetry  leads  to  the  total  spar  mass  matrix,  Nq^  x  Nq„ 
in  dimension,  given  by 
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tot 


(2.23) 


where  the  spar  cap  area  is  given  by  Eq.  1.3.  Using  Eq.  2.21  for  x,  each  element  of 
becomes  a  lineai'  combination  of  line  integrals  of  the  form 


Igp  {m,n) 


1 


jy=^yL 


cosAJ3'=% 


(51y  +  *S2) 


(2.24) 


Appendix  C  discusses  how  an  integral  table  for  different  m  and  n  combinations  may  be 
analytically  constructed  for  each  spar  line  integral. 

The  contribution  to  the  potential  or  strain  energy  U  of  an  element  of  length  ^^rl  of  the 


jsth  spar  is 


dU  =  -E  A  (ri) 


(2.25) 


where  E^^is  the  longitudinal  modulus  of  elasticity,  the  cap  area  is  a  function  of  Tj,  and 
;z(ri)  is  the  wing  depth  distribution  along  the  spar  line  (LIV90).  It  is  necessary  to  replace  T) 
dependence  with  y  dependence  in  this  expression.  Referring  to  Figure  A.3  and  Eq.  2.20 
the  transformation  for  the  bending  strain  is  given  by 

(^yR-^ytY 

The  cap  area  may  be  expressed  as  a  linear  function  of  y  as  shown  in  Eq.  1.3,  and  the  depth 
distribution  may  be  expressed  as  a  function  of  y  as  shown  in  Eq.  A.74  using  the  spar  line 
equation  of  Eq.  2.21.  Substituting  Eqs.  2.20  and  2.26  into  Eq.  2.25  gives 


(2.27) 


Using  Eq.  2.6  for  wix,y,t)  and  integrating  over  the  length  of  the  spar  gives 


U.  =  ^  ( cos  A)  ^Ej  ?  (y) 


hiy) 


{qj  («w’yy>  {qJ^V 


(2.28) 
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where  in  {ci„}  all  x  terms  are  replaced  by  Eq.  2.21  before  differentiation  with  respect  to  y. 
Summing  the  potential  energy  contributions  of  the  Nsy  spars  on  the  upper  wing  surface 
and  multiplying  by  2  because  of  wing  symmetry  leads  to  the  total  spar  stiffness  matrix, 
Nq^  X  Nq,,  in  dimension,  given  by 


Ns„ 


js=l 


(2.29) 


Each  element  of  becomes  a  linear  combination  of  line  integrals  of  the  form  Isp(m,n) 
as  defined  in  Eq.  2.24. 


2.5  Rib  Cap  Contribution  to  Mass  and  Stiffness  for  a  CRT  Zone 


The  contribution  to  the  kinetic  energy  T  of  an  element  of  length  dx  of  theyrth  rib  is 


(2.30) 


where  p^>  is  the  constant  material  density  of  the  rib,  and  the  cap  area  is  a  linear  func¬ 
tion  of  X  (LIV90).  Using  Eq.  2.6  for  w(x,y,t)  and  integrating  over  the  length  of  the  rib 
gives 


Tjr  =  .{«»>■  {aj  {aJ^UJdx 


(2.31) 


where  the  limits  of  integration  rxp  and  rx^  are  shown  in  Figure  A.4.  Summing  the  kinetic 
energy  contributions  of  the  Nry  ribs  on  the  upper  wing  surface  and  multiplying  by  2 
because  of  wing  symmetry  leads  to  the  total  rib  mass  matrix,  Nq„  x  Nq^  in  dimension, 
given  by 


Nr„ 


7>=1 


(2.32) 


where  the  rib  cap  area  is  given  by  Eq.  1.4.  With  y  equal  to  a  constant  yp,B  along  the  length 
of  the  rib,  each  term  of  becomes  a  linear  combination  of  line  integrals  of  the  form 


(2.33) 


Appendix  C  discusses  how  an  integral  table  for  different  m  and  n  combinations  may  be 
analytically  constructed  for  each  rib  line  integral. 

The  contribution  to  the  potential  or  strain  energy  U  of  an  element  of  length  dx  of  the 
yrth  rib  is 


du.^  =  2  V 


h  (x) 


rb,. 


w,ldx 


(2.34) 


where  Ej^  is  the  longitudinal  modulus  of  elasticity,  the  cap  area  is  a  linear  function  of 
X,  and  h{x)  is  the  wing  depth  distribution  along  the  rib  line  (LIV90).  Using  Eq.  2.6  for 
w{x,y,t)  and  integrating  along  the  length  of  the  rib  gives 


=  U  r 

2  jrjx==r 


'A(X) 


h{xy 


rh:. 


(?w>  {Qjdx 


(2.35) 


Summing  the  stiffness  contribution  of  the  ribs  on  the  upper  wing  surface  and  multi¬ 
plying  by  2  because  of  wing  symmetry  leads  to  the  total  rib  .stiffness  matrix,  Nq„  x  Nq„  in 
dimension,  given  by 


Nr„ 


jr=\ 


(2.36) 


where  h{x)  is  given  by  Eq.  1.1  holding  y  constant,  and  the  rib  cap  area  is  given  in  Eq.1.4. 
With  y  equal  to  a  constant  y^/g  along  the  length  of  the  spar,  each  element  of  becomes 
a  linear  combination  of  line  integrals  of  the  form  I^g{m,n)  as  defined  in  Eq.  2.33. 


2.6  Concentrated  Mass  Contribution  for  a  CPT  Zone 

Concentrated  masses  having  a  magnitude  of  Mj,.  may  be  designated  at  any  point 
{Xjc,yjc)  on  the  wing  structure.  The  contribution  to  the  kinetic  energy  T  of  theycth  concen¬ 
trated  mass  is 
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(LIV90).  Using  Eq.  2.6  for  w(Xj„yjj)  and  summing  over  the  Nc  masses  of  a  zone  leads  to 
the  total  concentrated  mass  matrix,  Nq„  x  Ng„  in  dimension,  given  by 


Nc 


jc=l 


(2.38) 


2.7  Load  Contributions  for  a  CPT  Zone 


The  generalized  load  vector  {P}  defined  in  Chapter  1  may  contain  contributions  from 
distributed  pressure  loads  and  concentrated  force  loads  both  acting  in  the  vertical  direc¬ 
tion.  A  distributed  pressure  load  over  a  trapezoidal  panel  area  may  be  defined  by  a  poly¬ 
nomial  load  function  with  terms  given  by 

♦  =  (2.39) 


The  contribution  to  the  work  g  of  a  load  distributed  over  an  area  dx  by  dy  and  working 
through  the  wing  deflection  is  given  by 


dQ  =  (x,  y)  w  (x,  y,  t)  dxdy 


(2.40) 


Using  Eq.  2.6  for  w{x,y,t)  and  integrating  over  the  panel  area  leads  to  the  generalized  load 
vector  contribution  of  the  distributed  load  given  by 


7V(t) 


j=\  yx 


(2.41) 


The  vector  contains  Nq„  terms,  and  its  ith  term  may  be  expressed  as  a  linear  combi¬ 

nation  of  area  integrals  of  the  form  lT^{m,n)  as  defined  in  Eq.  2.11,  where 


m  =  m^-  +  mq- 
n  =  n^.  +  nq. 


(2.42) 

(2.43) 


For  a  concentrated  force  load  of  magnitude  P  acting  at  the  point  {Xj^,yjc),  the  load  vector 
contribution  is  simply 
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{PJ  =P{aJ 


(2.44) 


where  {aj  is  evaluated  at  For  a  single  zone,  the  total  generalized  load  vector 

{P},„,  is  obtained  by  adding  the  load  contributions  from  all  distributed  and  concentrated 
loads  acting  on  the  zone. 


2.8  Displacement  and  Stress  Output 


The  CPT  zone  generalized  displacements  resulting  from  the  solution  of  Eq.  1.12  may 
be  used  to  calculate  the  static  displacement  and  stress  at  any  point  on  the  zone.  Choosing 
a  point  (x„,„yo,a)  from  an  output  grid,  the  vertical  deflection  w  at  that  point  may  be  deter¬ 
mined  simply  by  using  Eq.  2.6  where  {a„,}  is  evaluated  at  the  point.  Using  Eqs.  1.5,  2.3- 
2.5,  and  2. 17,  stresses  at  a  point  on  the  jlih  skin  layer  may  be  determined  from 


a 

XX 

XX 

^yy 

^  =  [0]/ 

£  [ 

ji 

^  ^xy  J 

=  -z[Q]ji[W]  {qj 


(2.45) 


where  [W]  is  evaluated  at  the  point.  Using  Eqs.  2.6  and  2.26,  the  uniaxial  stress  at  a  point 
(XoM^yoxt)  on  the  jsth  spar  may  be  determined  from 


T 

^js  =  =  -^^;.cos2A{a^,yy}  {qj 


(2.46) 


where  {a^,yy}  is  evaluated  at  the  point  (again,  the  x  powers  of  {a„}  are  replaced  by  Eq. 
2.21  before  the  differentiation  with  respect  to  y).  Similarly,  the  uniaxial  stress  at  a  point 
(Xout^yout)  on  the  jrth  rib  may  be  determined  from 


T 

^jr  ~  ~^^jr^^xx  ~  ~^^jr^^w’xx^ 


(2.47) 


where  is  evaluated  at  the  point. 

The  generalized  displacements  for  a  particular  natural  frequency  resulting  from  the 
solution  of  Eq.  1 . 14  may  be  used  to  calculate  the  mode  shape  displacement  at  any  point  on 
the  zone.  As  with  static  deflection,  the  modal  deflection  w  at  a  point  (x^„„yom)  n^ay  be 
determined  by  using  Eq.  2.6  where  {a„,}  is  evaluated  at  the  point. 


Chapter  3: 

Equivalent  Plate  Formulation  Using 
First  Order  Shear  Deformation  Plate  Theory 


3.1  Overview 

This  chapter  contains  the  First  Order  Shear  Deformation  Plate  Theoiy  based  formula¬ 
tion  of  the  mass,  stiffness,  and  generalized  load  contributions  needed  to  solve  for  the  gen¬ 
eralized  displacements  in  Eq.  1.11.  The  chapter  begins  with  a  discussion  of  the 
displacement  and  strain  representations  used  in  this  approach.  Following  this  is  the  deri¬ 
vation  of  the  wing  skin  contributions  to  the  mass  and  stiffness  of  the  system.  Next  dis¬ 
cussed  are  the  derivations  of  the  spar  cap  and  rib  cap  contributions  to  mass  and  stiffness. 
This  is  followed  by  the  derivations  of  the  spar  web  and  rib  web  contributions  to  mass  and 
stiffness.  Also  included  is  a  discussion  of  the  contribution  of  concentrated  masses  to  the 
mass  of  the  system.  Following  this  is  the  formulation  of  the  load  vector  from  concen¬ 
trated  force  load  and  distributed  pressure  load  conditions.  Finally  presented  are  the  meth¬ 
ods  used  to  calculate  displacement  and  stress  output  from  the  generalized  displacement 
solution.  Appendix  B  contains  a  more  detailed  derivation  of  mass  and  stiffness  contribu¬ 
tions  from  wing  structural  components. 

3.2  FSDPT  Displacements  and  Strains 

First  Order  Shear  Deformation  Plate  Theory  (FSDPT)  differs  from  Classical  Plate 
Theory  in  that  the  Kirchoff-Love  assumption  is  not  employed.  Rather,  it  is  assumed  that 
plane  sections  normal  to  the  plate’s  raid-surface  remain  plane  but  not  necessarily  normal 
to  that  surface  after  deformation.  Hence,  the  shear  strains  and  may  not  be  ignored  as 
they  are  in  CPT  (RE84,LIV194).  It  is  assumed  that  the  out  of  plane  displacements  are 
small,  and  there  may  be  in-plane  displacement  of  points  at  the  plate’s  mid-surface  since 
the  wing  is  not  assumed  to  be  symmetric  with  respect  to  the  x-y  plane.  As  with  CPT,  the 
transverse  normal  strain  is  negligible  since  the  transverse  deflection  does  not  vary 
through  the  thickness  of  the  plate,  thus  resulting  in  the  transverse  normal  displacement  w 
being  a  function  of  x  and  y  only. 
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Using  FSDPT  the  displacements  in  the  x,  and  z  directions  respectively  are  approxi¬ 
mated  by 


u  (x,  y,  z)  =  Uq  (x,  y)  +  (x,  y) 

V  (jc,  y,  z)  =  Vq  (jc,  y)  -i-  z^y  (x,  y) 
w(x,y,z)  =  WQ(A:,y) 


(3.1) 

(3.2) 

(3.3) 


where  Uq,  Vq,  Wq  are  the  jc,y,z  displacement  components,  respectively,  of  a  point  along  the 
reference  mid-surface,  and  and  X]/,,  are  rotations  of  a  line  element,  originally  perpendic¬ 
ular  to  the  mid-surface  plane,  about  the  y  and  x  axes  respectively  (RE84,LrV194).  The 
associated  strains  are  given  by 


du 


dx  dx  ^dx 


(3.4) 


dy  dy  dy 


(3.5) 


= 


9y'''8x  3y  dx  v3y  dx  j 


(3.6) 


dWr 


(3.7) 


dWr 

V  = 


(3.8) 


(LIV194,RE84).  These  displacement  and  strain  relationships  form  the  basis  for  the  mass 
and  stiffness  matrix  formulations. 

To  use  the  Ritz  method  described  in  Chapter  1,  it  is  necessary  to  approximate  each  of 
the  five  x,y,t  dependent  deformation  fields  by  polynomial  series  given  in  vector  form  as 


where 
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T 

UQ{x,y,t)  =  {aj(x,y)}  {q^(t)} 

(3.9) 

T 

VQ{x,y,t)  =  {a^ix,y)}  {q2it)} 

(3.10) 

T 

'^xix,y,t)  =  {a3(x:,y)}  {q^it)} 

(3.11) 

T 

^y(x,y,t)  =  {a^(x,y)}  {q^it)} 

(3.12) 

T 

WQix,y,t)  =  {a5(x,y)}  {^5(?)} 

(3.13) 

{aj(x,y)}  =  y  ,...} 

(3.14) 

and  similar  expressions  are  used  for  {a^y,  {<33}^,  and  {a5}^(LrV194).  The  column 
vectors  {^2}.  {?4}’  and  {^5}  contain  the  polynomial  coefficients  which  are  the 

generalized  displacements  in  the  Ritz  formulation.  These  polynomials  have  Nq^,  Nq2,  Nq^, 
Nq^,  and  Nq^  terms  respectively,  and  the  x  and  y  powers  are  predetermined. 

When  the  vectors  {^i},  [q2],  [qz],  {(Ia)-,  and  {^5}  are  combined  into  one  column  vec¬ 
tor  of  the  form 


'J*  7^ 

f  ^2’ ^3’ ^4’ ^5 ^ 


(3.15) 


the  M,v,w  components  of  deformation  may  be  written  as 


V 

w 


=  { [>5o(jc,y)] +z[*5i  (•^,y)]}  {^(0} 


(3.16) 


.  .  .  .  ,  j-  1  r  mS0(r,i)  nS0(r,i) 

where  Sq  and  Si  are  matnces  contammg  polynomial  terms  of  the  form  x  y 

and  respectively.  Partitioned  into  sub  vectors  of  dimension  1  x  Nq^ 

(n=l,2,...,5),  the  matrices  are  given  as 


a;  0  0  0  0 

[5o(x,y)]  =  0  a2  0  0  0 

0  0  0  0 


(3.17) 
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[^1  Uy)] 


0  0  03  0  0 

0  0  0  aj  0 
0  0  0  0  0 


where  both  are  3  x  in  dimension  with 

^^tot  =  +  ^^2  +  ^^3  + 


(3.18) 


(3.19) 


based  on  the  length  of  the  subvectors.  Note  that  different  orders  of  approximation  polyno¬ 
mials  are  allowed  for  different  displacement  fields. 


3.3  Skin  Contribution  to  Mass  and  Stiffness  for  a  FSDPT  Zone 


The  contribution  to  the  kinetic  energy  T  of  an  infinitesimal  skin  element  dx  by  dy  of 
thej/th  skin  layer  of  a  single  trapezoidal  panel  is 


T 

li 

li 

<  “ 

>  < 

~  > 

V 

V 

.  w  . 

.  w  > 

dxdy 


(3.20) 


where  p^,  is  the  constant  material  density  of  the  skin  layer,  and  tji{x,y)  is  the  thickness  of  the 
skin  layer  (LIV194).  Using  Eq.  3.16  for  the  displacements  and  integrating  over  the  whole 
trapezoidal  panel  area  (placed  at  a  distance  z  from  the  reference  plane)  gives 


(3.21) 


Summing  the  kinetic  energy  contributions  of  the  Mtot  skin  layers  in  the  panel  and  sum¬ 
ming  over  the  Np  panels  in  the  zone  leads  to  the  total  skin  mass  matrix,  Nqt^,  x  Nqi^,  in 
dimension,  given  by 

Np  NIjqj 

ip=\  jl=l  yx  (3.22) 


+  hji  {x,  y)  j  -I-  h^i  (x,  y)  (^6'[5  J  ]  JxJy 
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where  hji(x,y),  the  layer  depth  distribution  substituted  for  z,  and  tji{x,y)  are  polynomial 
functions  of  x  and  y  as  given  in  Eqs.  1.1  and  1.2  (LIV194).  Thus,  each  element  of 
is  a  linear  combination  of  area  integrals  of  the  form  Ij[i{m,n)  as  defined  in  Eq.  2.1 1. 

In  assessing  the  strain  energy  of  the  skin,  each  skin  layer  is  treated  as  a  plane  stress 
panel  for  which  the  only  strains  of  concern  are  and  Substituting  Eqs.  3.9-3.13 

into  Eqs.  3. 4-3. 6  allows  these  strains  to  be  written  as 


where 


Y 


xy 


{  [/?Q  (x, y)]  +z  [/?!  U, y)  ]  }  {q  (t)  } 


[/?Q  (x,  y)  ] 


[/?j  (x,y)] 


0  0  0  0 
0  al,y  0  0  0 
a[,y  al,^  0  0  0 

0  0  al,^  0  0 

0  0  0  al,y  0 

0  0  al,y  al,^  0_ 


(3.23) 


(3.24) 


(3.25) 


Partitioned  into  .subvectors  of  dimension  1  x  Nq„  (n=l,2,...,5),  these  matrices  are  3  x  Nqi^, 
in  dimension.  Now  the  contribution  to  the  potential  or  strain  energy  U  of  an  infinitesimal 
skin  element  dx  by  dy  of  the  jlth  skin  layer  of  a  single  trapezoidal  panel  is 


dUj,  = 


T 
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^vv 

XX 

XX 

'  £ 

£^v  ' 

yy 

yy 

Y  ^ 

*xv 

rdxdy 


(3.26) 


where  the  matrix  [Q]ji,  3  x  3  in  dimension,  is  the  jlth  layer’s  constitutive  matrix  referenced 
to  the  x-y  axes,  and  t,7(x,y)  is  the  thickness  of  the  skin  layer  (LIV194).  Using  Eq.  3.23  for 
the  .strains  and  integrating  over  the  trapezoidal  panel  area  gives 
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Vj,  =  {9)’'[<[eiy,«o+z(''o'2]/,) 

yx  (3.27) 

+2(  [  Q]  jiRq]  +  R]  [  G]  jiR  1 )  ]  {  ^ }  dxdy 

Summing  the  potential  energy  contributions  of  the  NIj-ot  skin  layers  in  the  panel  and  sum¬ 
ming  over  the  Np  panels  in  the  zone  leads  to  the  total  skin  stiffness  matrix,  Nq,^,  x  Nq,^,  in 
dimension,  given  by 

Np  ^^tot  (  t  ^ 

ip=\  y7=i  yx  (3.28) 

*h^,  ix,  y)  (/ibel,A)  +  '■),  U-  >■)  («[  muRijyxdy 

where  the  depth  distribution  hji{x,y)  in  Eq.  1.1  has  been  substituted  for  z,  and  tji(x,y)  is 
given  in  Eq.  1.2.  Each  element  of  is  a  linear  combination  of  area  integrals  of  the 

fonn 


3.4  Spar  Cap  Contribution  to  Mass  and  Stiffness  for  a  FSDPT  Zone 

The  contribution  to  the  kinetic  energy  T  of  an  element  of  length  dx[  of  theyyth  spar  is 


dT,,  = 
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' ,  ^ 
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.  vv  ^ 

dr\ 


(3.29) 


where  T|  is  the  coordinate  along  the  spar  axis  rotated  from  the  y  axis  by  an  angle  A,  is 
the  constant  material  density  of  the  spar,  and  the  cap  area  A  is  expressed  as  a  function 

R  js 

of  1)  (LIV194).  Referring  to  Figure  A.3,  all  T]  dependence  of  this  equation  may  be 
changed  to  y  dependence  using  Eq.  2.20.  The  spar  cap  area  may  be  expressed  as  a  linear 
function  of  y,  and  the  variable  x  may  be  expressed  in  terms  of  y  using  the  spar  line  equa¬ 
tion  from  Eq.  2.21.  Using  Eqs.  2.20  and  3.16  and  integrating  over  the  length  of  the  spar 
gives  p  .  . 

(3.30) 


2cosAJ>'=''>'z.  *-0  0  V  0  1 


-rz  sX  -rz"  5;5j]{^}(iy 
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Summing  the  kinetic  energy  contributions  of  the  NSfoT  spars  in  a  zone  leads  to  the  total 
spar  mass  matrix,  Nq,^,  x  Nq,„,  in  dimension,  given  by 


NSt 


=  2 


js=l 


o  +  hj,ix,y) 


<5, 


(3.31) 


where  the  spar  depth  distribution  hjXx,y)  has  been  substituted  for  z  using  Eq.  1.1,  and  the 
spar  cap  area  is  given  in  Eq.  1.3  (LIV194).  When  Eq.  2.21  is  used  for  x,  each  element  of 
VAIsp\tot  becomes  a  linear  combination  of  spar  line  integrals  of  the  form  Isp{m,n)  as  defined 
in  Eq.  2.24. 

The  contribution  to  the  potential  or  sti'ain  energy  (/  of  an  element  of  length  dy\  of  the 
jsth  spai‘  is 


‘‘Vjs  =  i  V 


(3.32) 


where  Ej,  is  the  longitudinal  modulus  of  elasticity,  the  cap  area  is  a  function  of  the 
spar  line  coordinate  r\,  and  is  the  normal  strain  along  the  spar  axis  (LIV194).  The 

spar  cap  area  may  be  expressed  as  a  linear  function  of  y.  Referring  to  Figure  A.3  and 
using  standard  tensor  transformation  rules,  the  normal  strain  in  the  r|  direction  may  be 
written  in  terms  of  strains  in  the  x  and  y  directions  by 


'xy 


(3.33) 


where 


{TR}^  =  {  sin^A,  cos^A,  sinAcosA}  =  {^^A,  c^A,  ^AcA} 


(3.34) 


Let  us  define 

[Qsp]  =  {TR}  {TR}"' 


(3.35) 


Substituting  Eq.  2.20  for  dr\,  using  Eqs.  3.23  and  3.33  for  the  strain,  and  integrating  over 
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the  length  of  the  spar  gives 


2cosAJy=^yL 


+z('f[  IQ„]  «o)  [G,,l  K,)]Mdy 


(3.36) 


Summing  the  potential  energy  contributions  of  the  spars  in  the  zone  leads  to  the  total 
spar  stiffness  matrix,  Nq,„,  x  in  dimension,  given  by 


P=l  (3.37) 

+'■„  7’)  ( [e„l  )  +  >>1  :p)  (  «[  [ Qj  ''l )  ]  ‘'7’ 


where  the  spar  depth  distiibution  hjXx,y)  has  been  substituted  for  z  using  Eq.  1.1,  and  the 
spar  cap  area  is  given  in  Eq.  1.3  (LIV194).  When  Eq.  2.21  is  used  for  x,  each  element  of 
[Kpliot  becomes  a  linear  combination  of  spar  line  integrals  of  the  form  Isp(m^). 


3.5  Rib  Cap  Contribution  to  Mass  and  Stiffness  for  a  FSDPT  Zone 

The  contribution  to  the  kinetic  energy  T  of  an  element  of  length  dx  of  theyrth  rib  is 


dT.^  = 
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(3.38) 


where  pj^  is  the  constant  material  density  of  the  rib,  and  the  cap  area  is  a  linear  func¬ 
tion  of  X.  Using  Eq.  3.16  for  the  displacements  and  integrating  over  the  length  of  the  rib 
gives 


T 


jr 


(3.39) 


where  the  limits  of  integration  rxp  and  rx^  are  shown  in  Figure  A.4.  Summing  the  kinetic 
energy  contributions  of  the  Nrj-oT  ribs  in  the  zone  leads  to  the  total  rib  mass  matrix,  Nqi^i  x 
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in  dimension,  given  by 


tot 


(3.40) 


where  the  rib  depth  distribution  hjXx,y)  has  been  substituted  for  z,  and  the  rib  cap  area  is 
given  by  Eq.  1 .4.  With  y  equal  to  a  constant  y^B  along  the  length  of  the  spar,  each  element 
of  [Mrhltot  becomes  a  linear  combination  of  line  integrals  of  the  form  InsifnA)  as  defined  in 
Eq.  2.33. 

The  contribution  to  the  potential  or  strain  energy  f/  of  an  element  of  length  dx  of  the 
jnh  rib  is 


dUjr  =  ^E.A  {X) 


rh/xx^^ 


(3.41) 


where  Ej^  is  the  longitudinal  modulus  of  elasticity,  and  the  cap  area  is  a  linear  func¬ 
tion  of  X.  Focusing  on  a  subset  of  the  displacement  fields  (using  Eqs.  B.4,  B.9,  and  B.ll) 
the  uniaxial  strain  along  a  rib  cap  is  expressed  by 

^XX  ^  ^  +z[Y^  (A:,y)]  }  {qit)} 

(3.42) 


where 


{q}^  =  {ql,ql} 


(3.43) 


[Eo(x,y)] 


{0} 


(3.44) 


[Tj  (x,y)] 


{0)  {4.,} 


(3.45) 


Both  matrices,  1  x  {Nq^+Nq-X  in  dimension,  are  partitioned  into  2  subvectors  having  Nq^ 
and  Nq-^  terms  respectively.  Using  Eq.  3.42  for  the  strain  and  integrating  over  the  length  of 
the  rib  gives 
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Summing  the  potential  energy  contributions  of  the  Nr-f^f  ribs  in  the  zone  leads  to  the  total 
rib  mass  matrix,  (Nqi+Nq^)  x  (Nqi+Nq^)  in  dimension  (and  submatrix  entries  correspond¬ 
ing  to  {q^}  and  {(73}),  given  by 


7>=i  (3.47) 

+v  >■)  ( >'o) + '>>  ( I'hi )  ] 

where  the  rib  depth  distribution  hjXx,y)  has  been  substituted  for  z,  and  the  rib  cap  area  is 
given  by  Eq.  1.4.  With  y  equal  to  a  constant  y,(,B  along  the  length  of  the  spar,  each  element 
of  {K^X\tot  becomes  a  linear  combination  of  line  integrals  of  the  form 


3.6  Spar  Web  Contribution  to  Mass  and  Stiffness  for  a  FSDPT  Zone 


A  spar  web  is  positioned  in  the  vertical  plane  between  parallel  spar  caps  on  the  upper 
and  lower  wing  surfaces.  This  plane  is  defined  to  be  the  r\-z  plane  where  r|  is  the  coordi¬ 
nate  along  the  corresponding  parallel  upper  and  lower  spar  axes  rotated  from  the  y  axis  by 
an  angle  A.  See  Figure  A.  3  for  spar  line  geometry.  The  contribution  to  the  kinetic  energy 
T  of  an  infinitesimal  element  dV[  by  dz  of  the  jlth  spar  web  layer  is 


dT. 
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dr[dz 


(3.48) 


where  py;  is  the  constant  material  density  of  the  web  layer,  and  tjiij\,z)  is  the  thickness  of 
the  layer  (LIV194).  The  layer  thickness  may  be  expressed  as  a  linear  function  of  y  only, 
and  the  variable  x  may  be  expres.sed  in  terms  of  y  using  the  spar  line  equation  from  Eq. 
2.21.  All  Ti  dependence  may  be  changed  to  y  dependence  using  Eq.  2.20.  Substimting  Eq. 
2.20  for  dr\,  using  Eq.3. 16  for  the  displacements,  and  integrating  over  the  area  of  the  web 
gives 
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2cosA« 


y=Sy^J"'L 


(3.49) 


T  \  if  T 

+d5i5j+z  5,5, 


]  {^}  dzdy 


where  the  limits  of  z  integration  are  the  depth  distributions  of  the  lower  and  upper  spar 
caps  given  by  Eqs.  B.133  and  B.134.  Summing  the  kinetic  energy  contributions  of  the 
spar  web  layers  of  each  of  the  Nsw  spar  webs  in  the  zone  leads  to  the  total  spar  web 
mass  matrix,  Nq,„,  x  in  dimension,  given  by 


[M 


1  _  r  jl  p-5yy(  f/ly 
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is=l  jl=l 


o  +  4Vi 
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(3.50) 


where  tjiiy)  is  given  in  Eq.  B.130.  The  z  integration  may  be  performed  analytically  using 
Eq.  B.  137  thus  leaving  only  a  spar  line  integral  to  be  evaluated.  Therefore,  when  Eq.  2.21 
is  used  for  x,  each  element  of  becomes  a  linear  combination  of  line  integrals  of  the 

form  Isp{m,n)  as  defined  in  Eq.2.24. 

In  assessing  the  strain  energy  of  a  spar  web,  each  layer  is  treated  as  a  plane  stress 
panel  where  the  only  strains  of  importance  are  and  From  the  assumptions 

of  FSDPT,  may  be  neglected.  In  Eq.  3.33,  the  normal  strain  has  been  defined  in 
terms  of  strains  in  the  jc-y  plane.  Using  standard  tensor  transformation  rules,  the  sheai- 
strain  y  may  also  be  defined  in  terms  of  strains  in  the  x-y  plane  by 

*  1 4 


=  {sinA,  cosA}  {sA,  cA}  |  [ 

^  ^  r,  J  I  ly,  >  (3.51) 

where  A  is  the  angle  of  rotation  of  the  T)  axis  from  the  y  axis  (LIV194).  Combining  Eqs. 
3.33  and  3.51  let  us  define 


where 


{E}  -  {^xx^^yy^'^xy’ 


(3.52) 


(3.53) 
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[7^^251  = 


(j'^A)  (c^A)  (i'AcA)  0  0 

0  0  0  (^A)  (cA)_ 


Using  Eqs.  3.4-3.15  to  expand  Eq.  3.53  gives 
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(3.54) 


(3.55) 


(3.56) 


(3.57) 


Paititioned  into  subvectors  of  dimension  1  x  Nq^  («=1,2,...,5),  both  matrices  are  5  x  Nq,^, 
in  dimension.  Now  the  contribution  to  the  potential  or  strain  energy  U  of  an  infinitesimal 
element  dr]  by  dz  of  the  jlth  spar  web  layer  is 


^r\z 


[S]„ 


■"mi 


'TIZ 


d'x\dz 


(3.58) 


where  [Sly/  is  the yVth  layer’s  constitutive  matiix  referenced  to  the  ri-z  axes  as  discussed 
in  Section  1.5,  and  r^7(Ti,z)  is  the  thickness  of  the  layer  (LIV194).  The  layer  thickness  may 
be  expressed  as  a  lineai-  function  of  y  only,  and  the  variable  x  may  be  expressed  in  terms  of 
y  using  the  spar  line  equation  from  Eq.  2.21.  Let  us  now  define 
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(3.59) 


Substituting  Eq.  2.20  for  dr[,  combining  Eqs.  3.52,  3.55,  and  3.59  for  the  strains,  and  inte¬ 
grating  over  the  area  of  the  web  gives 


t/.-,  =  I-  ^ 


2  cos  A 


+Z 


{q}  dzdy 


(3.60) 


where  the  limits  of  z  integration  are  the  depth  distributions  of  the  lower  and  upper  spar 
caps  given  in  Eqs.  B.133  and  B.134.  Summing  the  potential  energy  contributions  of  the 
spar  web  layers  of  each  of  the  Nsw  spar  webs  in  the  zone  leads  to  the  total  spar  web 
stiffness  matrix,  x  in  dimension,  given  by 


Nsw^Kwb 

L 

is=i  jl~\ 


+z(  w[  [ Q^]  Wq)  +  z^[  w]  [Q^]  Wj ]  ]  dzdy 


(3.61) 


where  tjiiy)  is  given  in  Eq.  B.130.  The  z  integration  may  be  performed  analytically  using 
Eq.  B.  137  thus  leaving  only  a  spar  line  integral  to  be  evaluated.  Therefore,  each  element 
of  becomes  a  linear  combination  of  line  integrals  of  the  form  Isp(m,n). 


3.7  Rib  Web  Contribution  to  Mass  and  Stiffness  for  a  FSDPT  Zone 


A  rib  web  is  positioned  in  the  vertical  plane  between  parallel  rib  caps  on  the  upper 
and  lower  wing  surfaces.  This  plane  is  the  x-z  plane  located  at  y-ym-  The  contribution  to 
the  kinetic  energy  T  of  an  infinitesimal  element  dx  by  dz  of  they/th  rib  web  layer  is 
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(3.62) 


where  p^/  is  the  constant  material  density  of  the  web  layer,  and  tji{x,z)  is  the  thickness  of  the 
layer.  The  layer  thickness  may  be  expressed  as  a  linear  function  of  x  only.  Using  Eq.  3. 16 
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for  the  displacements  and  integrating  over  the  area  of  the  web  gives 


J  ]  {4}  dzdx 


(3.63) 


where  the  limits  of  z  integration  are  the  depth  distributions  of  the  upper  and  lower  rib  caps 
given  respectively  by  Eqs.  B.  133  and  B.  134.  Summing  the  kinetic  energy  contributions  of 
the  of  rib  web  layers  of  each  of  the  Nrw  rib  webs  in  the  zone  leads  to  the  total  rib  web 
mass  matrix,  x  Nq,„,  in  dimension,  given  by 


tot 
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1 1  Pri 

ir=l  jl=l 


■x=rx^  rhy  (y) 
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(3.64) 


where  tjiix)  is  given  in  Eq.  B.199.  The  z  integration  may  be  performed  analytically  using 
Eq.  B.137  thus  leaving  only  a  rib  line  integral  to  be  evaluated.  Therefore,  with  y  equal  to 
a  constant  a  along  the  rib  line,  each  element  of  becomes  a  linear  combination  of 

line  integrals  of  the  form  Ineiniqi)  as  defined  in  Eq.2.33. 

In  assessing  the  strain  energy  of  a  rib  web,  each  layer  is  treated  as  a  plane  stress  panel 
where  the  only  strains  of  importance  are  e^,  and  y^.  From  the  assumptions  of  FSDPT, 
may  be  neglected.  Using  Eqs.  3.4,  3.7,  3.9,  3.11,  and  3.13  let  us  define 


where 
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(3.65) 
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[Fj  (x,y)] 


0  al,^  0 
_0  0  0 


(3.68) 


Partitioned  into  subvectors  of  dimension  1  x  Nq„  (n=l,3,5),  both  matrices  Fq  and  Fj  ai'e 
2  X  (A^^,+A^^3+A^^5)  in  dimension.  Now  the  contribution  to  the  potential  or  strain  energy  U 
of  an  infinitesimal  element  dx  by  dz  of  the  jlth  rib  web  layer  is 


dUji  =  ^tji(x,z) 


dxdz 


(3.69) 


where  is  the  jlth  layer’s  constitutive  matrix  referenced  to  the  x-z  axes  as  discussed 

in  Section  1.5,  and  tji(x,z)  is  the  thickness  of  the  layer.  The  layer  thickness  may  be 
expressed  as  a  lineai'  function  of  x  only.  Using  Eq.  3.65  for  the  strains  and  integrating  over 
the  area  of  the  web  gives 


(3.70) 


+z(  f[ [ 0] jiF^ ]  +  z\F\[mjiFM{q}  dzdx 


where  the  limits  of  z  integration  are  the  depth  distributions  of  the  upper  and  lower  rib  caps 
given  in  Eqs.  B.133  and  B.134.  Summing  the  potential  energy  contributions  of  the  of 
rib  web  layers  of  each  of  the  Nrw  rib  webs  in  the  zone  leads  to  the  total  rib  web  stiff¬ 
ness  matrix,  {Nq^+Nq^i+Nqj)  x  {Nq^+Nq^+Nq^)  in  dimension,  given  by 
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ir=l  jl=\ 


hLiy)  F 


ix)  [Fl[mjiF,^z[F,[mjiF, 


^z[f\  +  z"(f[  VmjiF^^yzdx 


(3.71) 


where  tji{x)  is  given  in  Eq.  B.199.  The  z  integration  may  be  performed  analytically  using 
Eq.  B.137  thus  leaving  only  a  rib  line  integral  to  be  evaluated.  Therefore,  each  element  of 
becomes  a  linear  combination  of  line  integrals  of  the  form  /^5(m,n). 
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3.8  Concentrated  Mass  Contribution  for  a  FSDPT  Zone 

Concentrated  masses  having  a  magnitude  of  Mj^  may  be  designated  at  free  moving 
points  {Xj,,,yjc)  on  the  wing  structure.  Each  point  is  also  associated  with  some  depth  distri¬ 
bution.  The  contribution  to  the  kinetic  energy  T  of  theycth  point  mass  is 

T  ■ 

u 

<  ” 

V 

.  vv 

Using  Eq.  3.16  for  the  displacements  and  summing  over  the  Nc  masses  of  a  zone  leads  to 
the  total  point  mass  matrix,  Nq,„,  x  Nq,^t  in  dimension,  given  by 


(3.72) 


=  X%[  Vo 


z\S,S,Uz 


(3.73) 


where  5o  and  5,  are  evaluated  at  ixjc,yjc)-  The  z  coordinate  of  a  point  is  determined  from  its 
depth  distribution. 


3.9  Load  Contributions  for  a  FSDPT  Zone 

As  discussed  in  Section  2.7,  the  generalized  load  vector  {P}  may  contain  contribu¬ 
tions  from  distributed  pressure  loads  and  concentrated  force  loads.  A  distributed  load 
over  a  trapezoidal  panel  ai'ea  was  defined  in  Eq.  2.39.  For  a  FSDPT  zone  the  pressure  may 
have  components  in  the  x  and  y  directions  as  well  as  the  z  direction.  If  the  component  dis¬ 
tributions  in  the  x,  y,  and  z  directions  are  (j)^,  (j)^,  (j)^  respectively,  then  the  3  directional  con¬ 
tributions  to  the  potential  energy  Qof  a  load  distributed  over  an  area  dx  by  dy  of  a  panel 
and  working  through  the  wing  deflections  are  given  by 


dQ^  =  {x,  y)  u  (x,  y,  z,  t)  dxdy  =  (x,  y)  [mq  (x,  y,  t)  +  z^^f^  (x,  y,  t)  ]  dxdy 
dQy  =  ^y{x,y)v{x,y,z,t)dxdy  =  (t)j,(x,y)  [vQ(x,y,  t)  +zyfy{x,y,t)]dxdy 
dQ^  =  (^^(x,y)w(x,y,z,t)dxdy  =  (|)^(x,  y)  [wQix,y,  t)]  dxdy 


(3.74) 

(3.75) 

(3.76) 


Using  3.9-3.13  for  the  deformation  fields,  we  get  an  expression  for  the  virtual  work  in  the 
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form 


xy 

%z{a^}  {6^4}  +(l)Ja5}  {5^5}  yxdy 
Let  the  polynomial  loads  now  be  defined  as 


(3.77) 


tmifj  ni^j 


=  I V  ’V 

j=  1 


m^j  niS)j 


=  I V  V 

;=i 


mfyj  n<S)j 


=  I V  V 

;■  =  1 

The  generalized  load  vector  is  now  partitioned  as 


T  T  T  T  r 
{P}'  =  {P^,P^,PyP^,P^} 


where 


j=\  yx 


{P2}  = 

7=1 


{  ^3  }  =  ^  X  i  {  a3  (JC,  y)  }  dxdy 

7=1  3'^ 


N(?, 


^^4}  =  ^I%jKV^a4(7C,y)}./x^y 

7=1 


7=1 


(3.78) 


(3.79) 


(3.80) 


(3.81) 


(3.82) 


(3.83) 


(3.84) 


(3.85) 


(3.86) 
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The  depth  distribution  of  the  loaded  panel  is  substituted  for  z  in  terms  where  it  appears, 
depending  on  where  the  in-plane  forces  are  applied.  Each  term  of  the  vector  may  be 
expressed  as  a  linear  combination  of  area  integrals  of  the  form  as  defined  in  Eq. 

2. 1 1.  For  a  concentrated  force  load  with  components  P^,  Py,  P^  acting  at  point  the 

load  vector  contribution  is  simply 


iPc}  = 


^ 1 } 
Py{a2} 
PxZ{a^} 
PyZ{a^} 
Pzi^5^ 


(3.87) 


where  all  polynomial  terms  are  evaluated  at  and  z  is  evaluated  from  the  points 

depth  distribution.  For  a  single  zone,  the  total  load  vector  {P},^,  is  obtained  by  adding  the 
load  contributions  from  all  distributed  and  concentrated  loads  acting  on  the  zone. 


3.10  Displacement  and  Stress  Output 

The  FSDPT  zone  generalized  displacements  resulting  from  the  solution  of  Eq.  1.12 
may  be  used  to  calculate  the  static  displacement  and  stress  at  a  selected  point  in  the  zone. 
Choosing  a  point  (x„„t,yout)  from  an  output  grid  and  determining  Zom  from  the  point’s  depth 
distribution,  the  displacements  u,  v,  and  w  at  that  point  may  be  determined  by  using  Eq. 
3.16  where  So  and  Sj  are  evaluated  at  the  point.  Using  Eqs.  1.5  and  3.23,  stresses  at  a  point 
(x„„i,yoJ  on  the  Jlth  skin  layer  may  be  determined  from 


yy 


Jty  ji 


=  [Q]ji  {  [Rq  U  y)  ]  +  z^^i  [/?!  (x,  y)  ]  }  {g} 


(3.88) 


where  Rq  and  are  evaluated  at  the  point,  and  z<,„,  is  evaluated  from  the  layer’s  depth  dis¬ 
tribution  at  the  point.  Using  Eqs.  3.23  and  3.33,  the  uniaxial  stress  at  a  point  (x:„,„,yoJ  on 
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the  jsth  spar  may  be  determined  from 

=  Ej.  { ™}  { ['*0  y)  1  +  [«1  (^.  y)  ] }  { 9} 

where  Rq  and  are  evaluated  at  the  point,  and  Zout  is  evaluated  from  the  spar’s  depth  dis¬ 
tribution  at  the  point.  Using  Eq.  3.42,  the  uniaxial  stress  at  a  point  (x^„„yout)  on  theyrth  rib 
may  be  determined  from 


(Sjr  =  {  [  Yq  (x,  y)  ]  +  U  y)  ]  }  {  9) 


(3.90) 


where  Yq  and  are  evaluated  at  the  point,  and  Zom  is  evaluated  from  the  rib’s  depth  distri¬ 
bution  at  the  point. 

The  generalized  displacements  for  a  particular  natural  frequency  resulting  form  the 
solution  of  Eq.  1.14  may  be  used  to  calculate  the  mode  shape  displacement  at  a  selected 
point  in  the  zone.  As  with  static  displacements,  the  modal  displacements  at  a  point 
y^ut)  may  be  determined  by  using  Eq.  3.16  where  So  and  Si  are  evaluated  at  the  point, 
and  Zout  is  evaluated  from  the  point’s  depth  distribution. 


Chapter  4: 

Zone  Connections  and 
Boundary  Conditions 


4.1  Overview 

This  chapter  presents  the  global  synthesis  of  a  multi-zone  wing  model.  The  chapter 
begins  by  discussing  how  the  stiffness,  mass,  and  load  contributions  from  multiple  CPT 
and  FSDFPT  zones  are  merged  into  the  global  math  model.  Next  presented  is  the  deriva¬ 
tion  of  stiffness  matrix  contributions  from  the  spring  connection  of  two  CPT  zones.  Fol¬ 
lowing  this  is  the  derivation  of  stiffness  matrix  contributions  from  the  spring  connection 
of  two  FSDPT  zones.  These  are  extended  with  the  derivation  of  the  stiffness  matrix  con¬ 
tributions  from  the  spring  connection  of  a  CPT  zone  and  a  FSDPT  zone.  Finally  discussed 
is  the  imposition  of  wing  root  boundary  conditions. 


4.2  Zone  Connections  in  Global  Context 

When  multiple  zones  are  used  to  characterize  a  wing,  compatibility  of  displacements 
at  the  zone  interfaces  must  be  enforced.  This  may  be  done  using  a  penalty  method  via 
computational  springs.  A  spring  stiffness  coefficient  represents  local  stiffness  at  an  attach 
point  common  to  two  zones.  A  large  stiffness  coefficient  represents  very  rigid  connection 
while  a  smaller  coefficient  may  represent  more  flexible  connection,  such  as  that  provided 
by  an  actuator.  Linear  springs  may  be  used  to  enforce  compatibility  of  linear  displace¬ 
ments,  and  rotational  springs  may  be  used  to  enforce  compatibility  of  slopes  or  angular 
displacements. 

In  a  multiple  zone  configuration,  the  stiffness  matrix,  mass  matrix,  and  load  vector 
contributions  of  each  zone  must  be  merged  into  their  global  counterparts, 
and  {P}giob-  These  correspond  to  the  vector  of  global  generalized  displacements  formed 
from  the  generalized  displacements  of  each  zone.  As  an  example,  a  wing  with  2  FSDPT 
zones  (zones  #1  and  #2)  and  2  CPT  zones  (zones  #3  and  #4)  will  have  a  vector  of  global 
generalized  displacements  given  by 
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T  T  T  T  T  T 

glob  ~  ^1,’ ^2,’ ^3,’ ^4,’ ^5, 


^2^’  ^3/ ^42’ 


r 

^w, 


^  1 


(4.1) 


(Refer  to  Chapters  2  and  3  for  the  notation  of  generalized  displacements  for  CPT  and 
FSDPT  zones).  The  additional  second  subscript  here  denotes  the  zone  number  out  of  a 
possible  Nz  zones.  The  vector  in  Eq.  4.1  is  of  length  Nqgi^i,  where 


^^gloh  = 


(4.2) 


Thus,  the  global  stiffness  and  mass  matrices  will  be  x  in  dimension,  and  the 
global  load  vector  will  be  of  length  Nq^i^^.  The  global  stiffness  matrix  may  be  partitioned 
into  submatrices  as  follows: 


[K 


glob^ 


^1 1 

^1  2 

IjZj 

3 

^14 

^1  5  1 

1 

2 

M^2 

K,  3 

^14 

^1,5: 

^2  1 
Zj  ij 

^2  2 

^2  3 

Zj 

^2  4 

ZjH^j 

^2  5  1 

Zj  Jj  1 

^2  1 

Zj  i2 

^22 

ZjZ2 

^2  3 

Zj32 

^2  4 

^2,5, 

^3  1 

Ojlj 

^3  2 

^3  3 

^3  4 

^3  5  1 

jj  1 

^3  1 
^2 

^3  2 

^33 

^3  4 
■’1^2 

^2,5, 

^4  1 

M 

^4  2 

^4  3 

^4  4 

^45 1 

^4  1 

^2 

^4  2 

^43 

^4  4 

^4,5. 

^5  1 

Jj  ij 

^5  2 

JjZj 

^5  3 

^5  4 

^55 1 

jj  j,  1 

^5  1 

■^1^2 

^5  2 
Jl'^2 

^5  3 

^5  4 

-’1^2 

^5.5. 

1 

■‘2M 

^12 

3 

^14 

^15  1 

^1  1 

•^2^2 

^1  2 

/f,  3 

^2^2 

^14 

^2^2 

^2  1 

Z^l  1 

^22 

^2  3 

^2  4 

^2  5  1 

^2-^1  1 

^2  1 

‘^2-‘-2 

^22 

^2^2 

^2  3 

^2'^2 

^2  4 

^2^2 

^3  1 

^3  2 

^33 

^3  4 
•^2^1 

^3  5  1 

^2-^1  I 

^3  1 

^2^2 

^3  2 
^2^2 

^3  3 
^2^2 

^34 

^2^2 

^4 1 

^4  2 

^4  3 

^4  4 

^4  5  1 
^2-^1  1 

^4  1 

^2^2 

^4  2 

^2‘^2 

^4  3 
^2^2 

^4  4 
^2^2 

^4.5^ 

^5  I 

^5  2 

^5  3 

^5  4 

•^2^1 

^5  5  1 

■^2'’l  1 

^5  1 

^2^2 

^5  2 

■^2^2 

^5  3 
^2^2 

^5  4 

•^2^2 

^Wjl, 


^M'42, 


^*>'342  ^*^352 


^W^22 


^2, .3 

^2..3 

^3..3 

^4..3 

^5..3 

1  ^23-^3! 

1  "^33-3! 

1  ^4,.3| 

1  ^53^3! 

^2’^4 

1  K 

^  W3 

1 

1 

^4>^3 

1 

(4.3) 


where  each  submatrix  ,  is  of  dimension  Nq^  x  Nq^  (r,.s=w,  1,2,3 ,4,5  and  /j=l,2,...,iVz). 

r^Sj  .  .  .  ‘ 

The  global  mass  matrix  will  be  partitioned  like  the  global  load  vector  will  be 

partitioned  like  {^}giob-  Contributions  to  stiffness  and  mass  from  each  zone’s  structural 
components  are  placed  in  the  corresponding  submatrices  on  the  diagonal  where,  for  a  sub¬ 
matrix  K .  ,  i  -j  -  zone  #.  There  are  no  mass  contributions  in  submati'ices  where  jVy . 

;  iSj 

However,  spring  connections  do  provide  connectivity  between  zones  and  therefore  pro¬ 
duce  contributions  to  stiffness  in  submatrices  where  i^j. 
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4.3  Spring  Connection  of  2  CPT  Zones 


Let  us  begin  by  defining  an  “attach  line”  which  corresponds  to  the  common  surface  of 
2  CPT  zones  to  be  connected.  Springs  act  at  selected  “attach  points”  along  this  line  to  pro¬ 
vide  local  stiffness.  The  left  and  right  endpoints  of  the  attach  line  are  given  by  the  coordi¬ 
nates  (atXi,atyi)  and  {atXi{,atyg).  Referring  to  the  attach  line  in  Figure  4.1,  let  us  define  the 
unit  vector  1  in  the  direction  of  the  attach  line  given  by 


1  =  cos  a/  +  sin  ay 


(4.4) 


where 


cos  a  = 


atXf^  -  atx^ 


sin  a  = 


jJ  (atXjf-atx^)^  +  {aty^  -  atjj) 

_ atyR  -  atyj^ _ 

JiatXj^  -  atXj)  ^  +  { aty^  -  atyj) " 


(4.5) 


(4.6) 


A  vector  of  rotation  having  positive  rotational  vector  components  in  the  x  and  y  directions 
may  be  given  by 


n  =  (QjU  {a/j 


(4.7) 


Figure  4. 1:  CPT  Attach  Line  Geometry 
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Using  the  right  hand  rule  for  a  CPT  zone,  corresponds  to  the  displacement  slope 
while  Qy  corresponds  to  -w,^.  Positive  rotation  about  i  may  then  be  given  as 

0  =  n»t  =  w,,cosa- w,^sina 

^  (4.8) 


(LIV190).  Recall  from  Section  1.3  that  each  zone  may  have  its  geometry  defined  in  an 
independent  axis  system.  Therefore,  the  attach  line  common  to  both  zones  may  be  defined 
with  a  unique  orientation  ttj  in  zone  1  and  a2  in  zone  2.  This  results  in 

Gj  =  Wp^costtj  -  Wp^sinrij 
02  =  W2,yCosa2-W2,^sina2 


(4.9) 

(4.10) 


where,  from  Eq.  2.6,  we  have 


T 


T 

^2  =  {a^}  {q^} 


(4.11) 

(4.12) 


Now  let  a  linear  spring  with  stiffness  connect  point  (xi,y,)  on  zone  1  with  point  (X2,y2) 
on  zone  2  by  providing  vertical  displacement  compatibility.  The  contribution  to  the  poten¬ 
tial  energy  U  of  the  spring  resisting  vertical  displacement  is 


where 


(4.13) 


A  =  Wj  -  W2 


(4.14) 


Substituting  Eqs.  4.14, 4.11,  and  4.12  into  Eq.  4.13  gives 


T  T  T  T 

T  T  T  T 


(4.15) 


This  leads  to  the  stiffness  submatrix  contributions  given  by 


=  K{a^}  {«w,> 


(4.16) 


where  {a  )  and  {a  )  are  evaluated  at  (A:i,yi)  and  (x2,y2) ^spectively.  These  submatri- 

*-  W|  -*  ^2 

ces  are  merged  into  the  global  stiffness  matrix  having  the  form 


^W2W, 


(4.17) 


for  each  linear  spring  connection  that  is  used. 

Also  let  a  rotational  spring  with  stiffness  connect  point  (a:i,yi)  on  zone  1  with  point 
(Xz,y2)  on  zone  2  by  providing  slope  compatibility.  The  contribution  to  the  potential  energy 
U  of  the  spring  resisting  angular  displacement  is 


(4.18) 


Substituting  Eqs.  4.9-4.12  into  Eq.  4.18  gives 


-  cosa^ sina^  J 

T  T\ 

— cosoCj sinoCj  +  sin  oCj 

T/T  T 

■{q^  }  l^cosa^cosa^ia^^  y}  {a^^y}  -  cosa^sina^la^^  y}  J 

T  T  \ 

-cosa^sinajla^^  J  +  sinttj  sina^  {a^  J  J  j{q^} 
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J.T  T  (4.19) 

-{«  }  cosa,cosa,{a^^J  {a^^J  -cosa2sinai  {a  }  {a  J 


-cosajSina2{a^^J  {a^  y}  +  sinajSina2  J  J  \{q^} 


T/  T  T 

+  {q\  cos2a2{a  T{a  >  -  cosa2sina2  {a  }  {a  } 


T  T  \ 

— COS OC^ since.,  ^2  I 


This  leads  to  the  stiffness  subraatrix  contributions  given  by 


r  T  T 

=  cos2aj{a^^  ^}  - cosa,sinai  J 

T  T- 

-cosaisinaj  {a^  J 

"  -^9|^(^osaiCOsa2{fl^^  3,}  {a^^^y}  -  co&a^&ina^  {a^^  y}  {a^^  J 

T  .  . 

— cos CC2 since j  { {<3^  y}  "t  sincej  since2  j,}  ^^W2,x^ 


(4.20) 


r  r  T 

=  -Ucosa^cosa^ia^^y}  {a^^y}  -  cosa^sina^ia^^y} 


+smaiSina2{a^^  J  J 


-cosajSina2 


=h  cos2a2{a^^,y}  -  cosa^sina^{a^^y}  {a^^  J 


-cosa2sina2  {%,3,}  +  sin^oei 


where  {a^  y}  and  {a^^  J  are  evaluated  at  (xiji),  and  {a^^  y}  and  {a^^  J  are  eval¬ 
uated  at  (X2,y2)  (LIV190).  These  submatrices  are  merged  into  the  global  stiffness  matrix 
having  the  form  shown  in  Eq.  4.17  for  each  rotational  spring  connection  that  is  used. 
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4.4  Spring  Connection  of  2  FSDPT  Zones 


Point  compatibility  of  relative  vertical  displacements  between  2  FSDPT  zones  may  be 
enforced  with  a  linear  spring  in  the  same  fashion  as  it  is  done  for  CPT  zones.  Let  a  linear 
spring  with  .stiffness  connect  point  (xi,yi)  on  FSDPT  zone  1  with  point  (X2,y2)  on  FSDPT 
zone  2.  U.sing  Eqs.  3.3  and  3.13  let  us  define 


T 

=  {a5,> 


T 

W'2  =  {%} 


(4.21) 

(4.22) 


These  relationships  for  Wj  and  W2  may  be  used  in  Eq.  4.13  to  express  the  contribution  to 
the  potential  energy  of  the  spring  as 


~  2^5 1^  ^^5,^ 

T  T  T  T  -i 

{<^52^  ^^52^  '^^52^  J 


(4.23) 


This  leads  to  stiffness  submatrix  contributions  given  by 


^^525,^  yt  ~  ^5  ^*^52^ 


(4.24) 


where  {a.}  and  {a.)  are  evaluated  at  (x„yi)  and  (:)r2»>'2)  respectively.  These  submatri- 

^  1  ^2 

ces  are  merged  into  the  global  stiffness  matrix  having  the  form 


51 


glob^ 


^1 1 
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22 
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32 
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1  ^ 

1 

b 

2^2 

>32 
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2^2 

^42 

>32 

^4  4 
^2^2 

^4 

1^2 

^5  1 

^5, 

:2, 

^5, 

.3, 

^5 

>4. 

^5 

2^1 

1  ^5 

2 

^2 

^5 

>22 

^52 

>32 

^5  4 

“*2^2 

^5 

(4.25) 


for  each  vertical  spring  connection  that  is  used. 

Linear  in-plane  displacement  compatibility  between  2  FSDPT  zones  must  also  be 
enforced.  Referring  to  Eqs.  3. 1-3.2,  the  in-plane  displacements  u  and  v  both  consist  of  a 
linear  translation  (mo»Vo)  and  a  rotation  (\)4,\Ky).  There  are  2  ways  to  enforce  point  compat¬ 
ibility  of  these  translational  and  rotational  components.  A  coupled  pair  of  lineai*  springs, 
both  at  a  distance  h  from  the  mid-surface,  may  be  used  to  enforce  compatibility  of  both 
components.  Alternately,  a  linear  spring  may  be  used  at  the  mid-surface  to  enforce  trans¬ 
lational  compatibility,  and  a  separate  rotational  spring  may  be  used  to  enforce  rotational  or 
slope  compatibility.  Both  methods  become  equivalent  when  the  stiffness  coefficient  of  the 
rotational  spring  is  chosen  to  be  the  stiffness  coefficient  of  the  linear  springs  divided  by  h  . 
Since  the  latter  method  has  no  complicating  depth  distribution  dependence,  it  is  the  one 
used  and  presented  in  this  work. 

Referring  to  Figure  4.2,  let  us  define  the  unit  vector  f,  perpendicular  to  the  attach  line, 
given  by 

r  =  sinai-cosaj 

(4.26) 


where  Eqs.  4.5  and  4.6  define  sina  and  cosa.  Let  us  also  define  a  vector  of  mid-surface 
in-plane  displacement  with  components  in  the  positive  x  and  y  directions  given  by 


X  =  V  + V 


(4.27) 


Lineal'  displacement  perpendicular  to  the  attach  line  at  the  mid-surface  is  then  given  by 
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Figure  4.2:  FSDPT  Attach  Line  Geometry 


^  •  f  =  Mpsina- VgCosa 


(4.28) 


Given  that  the  attach  line  may  be  oriented  differently  in  each  zone’s  axis  system,  Eqs.  3.9 
and  3.10  may  be  used  to  express  linear  displacements  for  zone  1  and  zone  2  respectively 
by 


T  T 


Cl  = 

sinaj{aj  } 

-cosaj{a2}  {^2,^ 

(4.29) 

T 

T 

C2  = 

sina2{aj^}  {q^} 

-  cosa.,  {<^9  }  {Qo  } 

Z  Z>2  ^2 

(4.30) 

Now  let  a  linear  spring  with  stiffness  ki  act  perpendicular  to  the  attach  line  and  con¬ 
nect  point  (x„yi)  on  zone  1  with  point  (X2,y2)  on  zone  2.  The  contribution  to  the  potential 
energy  U  of  the  spring  resisting  in-plane  displacement  is 


(4.31) 


Substituting  Eqs.  4.29  and  4.30  into  Eq.  4.3 1  gives 
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r  T  T  T  T 

sin2aj{(?i}  {a^}  {a^}  {5,  }  -  cosajsina  {a^}  {a^}  {q^} 


cosaj sin ttj  {(72^1  {^2,^ 


+  cos^a 


1  {^2,^ 


-sinajSina2{^i^}  {a^}  {a^}  {^j  J  +  cosa2sinaj  {^j_}  {a^}  {a^}  {q2} 


+ cos  a  j  sin  a2 


T  T  T  T 

{q2  }  {(^2}  ^<^\)  {^i^}-cosajCosa2{^2,^  Uh) 

(4.32) 


-sinttjSinaj +  costtj sina2 {^j^}  {^i^)  {^2,^ 

T  T  TT 

+cosa2sinaj  {(J2J  {^22^  {^j_} -cosajCosa2  {^2^^  ^^22^^21^ 

T  T  T  T 

&in^a2{q^}  («i^)  {«i^}  {^1^} -cosa2sina2{^i^}  {<^1^}  {^2^}  {^2^^ 

T  T  T  T  - 

-cosa2sina2  {?2j  ^^2^^  +  cos2a2  {^22^  ^^22^  ■>^^22^  ■^^^22^ 

This  leads  to  stiffness  submatrix  contributions  given  by 


^2,1,  ^^2,2.  ‘ 


j^sin^ttj  {aj  }  l^-cosajsinaj  {aj  }  {a2} 

l^-costtjsinaj  {a2  }  {Oj  ^cos2aj{a2}  {^2} 


^1  i  ^12 
^2,12  ^2,22 


-sinajsina2  {aj  }  ^cosa2sinaj  {aj  }  {<22^} 

cosajSina2 {^2^}  {‘^i2^’n  r-cosajCosa2  {^I2,)  ^*^22^ 


(4.33) 


^11  ^12 

^2-^1  -^2^1 

Kr.  t  Kj  r. 


-sinajSina2{aj^)  j^cosajSina2  {02  } 

cosa2sinaj  {(22  }  (oj  iH  r-cosajCOsa2 {<222^ 
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^11  ^12 

^2^2  ^2^2 

^0  1  ^7  7 

^2^2  ^2-^2 


T- 

T- 

sir?-a2{a^}  {a^} 

-cosa2sina2  {<^2.) 

n 

T- 

-cosa2sina2  {^2^}  i^i^} 

cos^a,  {^9  }  {flo  } 

^  ^2  ^2 

where  {a,  }  and  {a,  )  are  evaluated  at  (xi,)'i),  and  {a,  }  and  {a^}  are  evaluated  at 
(X2,y2)-  These  submatrices  are  merged  into  the  global  stiffness  matrix  having  the  form 
shown  in  Eq.  4.25  for  each  translational  spring  connection  that  is  used.  Note  that  if  a  lin¬ 
ear  spring  of  stiffness  k2  is  used  to  resist  translation  parallel  to  the  attach  line,  Eq.  4.33 
may  be  used  to  find  the  stiffness  contributions  with  sina  replaced  by  cosa,  cosa  replaced 
by  -sina,  and  ki  replaced  by  ^2- 

Enforcement  of  rotational  or  slope  compatibility  is  done  in  a  manner  similar  to  that 
for  translational  compatibility.  Recall  that  a  vector  of  rotations  was  defined  in  Eq.  4.7. 
Using  the  right  hand  rule  for  a  FSDPT  zone,  Q;,  corresponds  to  -\[t^  and  Q,  correspond  to 
Using  Eqs.  4.4  and  4.7,  positive  rotation  about  the  attach  line  may  be  given  as 

0  =  n*?  =  xir^sina-xi/,  cosa 

^  (4.34) 

Given  that  the  attach  line  may  be  oriented  differently  in  each  zone’s  axis  system,  Eqs.  3. 1 1 
and  3. 12  may  be  used  to  write  rotational  displacements  for  zone  1  and  zone  2  respectively 
given  by 


T  T 

ej  =  sinaj{a3}  {^3  }  -  cosa,  {<24  }  {^4} 
T  T 

02  =  sina2{<23}  -  cos,a^{a^} 


(4.35) 

(4.36) 


Now  let  a  rotational  spring  with  stiffness  k^,  act  in  resistance  to  rotation  about  the  attach 
line  and  connect  point  (a:i,yi)  on  zone  1  with  point  (a:2,>’2)  zone  2.  The  contribution  to 
the  potential  energy  U  of  the  spring  is 


t't,  =  ^*3(61-62)^ 


(4.37) 


Substituting  Eqs.  4.35  and  4.36  into  Eq.  4.37  gives 
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t^.3  =  2^3 


r  T  T  T  T 

{ag^}  {a3^}  {<73^}  -  costtjsinaj  {^3  }  {a^  }  {a^  }  {q^  } 

T  T  T  T 

-cosa^sinaj  {^4  }  {a^}  {a^}  {^3  }  +  cos^aj  {^4  }  {a^}  {a^}  {q^} 

T  T  T  T 

-smajsina2{^3  }  {<23  }  {a^  }  {^3  }  +  costtjsinaj  [q^]  {a^^}  {a^}  {q^} 


T  T  T  1 

+cosajSina2  {^4  }  {a4  }  {<33  }  {^3  } -cosajCosa2  {^4  }  {a^}  {a^}  {q^} 

(4.38) 


T  T 

}  {^3  }  {«4  }  {?4  J 

I  ^1  ^2  ^2 

T  T 

,}  {«4,> 


T  T  T  T 

-sinajSina2{93^}  {a^}  {a^}  {^3  }  +  cosaj  sina2  {<?3^}  {03^}  {a4}  {q^} 


^2  -^2 
T 


+cosa.,sina,  {g,  }  {a.}  {a^}  {^3  } -cosa.cosaj  {?4  }  {<34  }  {a4  }  {124} 

T  T  T  T 

sin^(X2{^73  ^  ^  ^^3  ^  cos CX2  sin  (X2  1  ^^3  ^  ^^4^^  ^^^2^ 


^2  -^2  ^^2  -"2 

r  T 


-cosa2sina2  {<?4  }  {04  }  {a.  }  {q.  }  +cos,^a2{q^  }  {a^  }  {a.}  {qA 

^  ^  *0  0.  0.  *7.  2  2  2  2 


2  -’2  ^2  ^2 

r  r 

?4,> 


^33 

'*'3, 4' 

=  h 

^4  4 

'^4,4, 

sin^aj  {a3  }  {03  }  j^-cc 

-cosajsinaj  {a4  }  {03^}  ^ 


cosajsinaj  {^3  }  {^4  } 


cos^aj  {^4  }  {^4  } 


■^3,42  _ 


^4  3  ^^4  4 


j^-sinajsina2  {^3  }  {a3^}^  j^cosa2sinaj  {^3  }  {a^} 
j^cosajSina2  {a4  }  {<23^)^  j^-cosa^cosa2  {^4  }  {<34^} 
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(4.39) 


■^3  3 

/^4  3  ^4  4 


=  kn 


j^-sinajsina2{a3j  j^cosajSina2 {ag  }  {a^} 

j^cosa2sinaj  }  {a3^}^  |^-cosajCOsa2 {a4^}  } 


K,  .  K, 


^^2^2  3,42 


‘■423,  "'‘4,42 


= 


|^sin^a2{a3  }  {<^3^  1  [  -cosa2sina2  {ag^}  {^4^}^ 

|^-cosa2sina2{a4^}  {Og^}^  |^cos^a2 {04^}  {<34^} 


where  {a^}  and  {a^}  are  evaluated  at  (Xi,y]),  and  {a^}  ^^42^  are  evaluated  at 

{X2,y2).  These  subraatrices  ai'e  merged  into  the  global  stiffness  matrix  having  the  form 
shown  in  Eq.  4.25  for  each  rotational  spring  connection  that  is  used.  Note  that  if  a  rota¬ 
tional  spring  of  stiffness  k^  is  used  to  resist  rotation  about  an  axis  perpendicular  to  the 
attach  line,  Eq.  4.39  may  be  used  to  find  the  stiffness  contributions  with  sina.  replaced  by 
-coscL,  cosa  replaced  by  sina,  and  k^,  replaced  by  k^. 


4.5  Spring  Connection  of  a  FSDPT  Zone  and  a  CPT  Zone 

Let  us  define  zone  1  to  be  a  FSDPT  zone  and  zone  2  to  be  a  CPT  zone.  Point  compat¬ 
ibility  of  relative  vertical  displacements  between  the  two  zones  may  be  enforced  as  fol¬ 
lows.  Let  a  linear  spring  of  stiffness  k^  connect  point  (x:i,yi)  on  zone  1  with  point  (X2,y2)  on 
zone  2  by  resisting  linear  vertical  displacement.  The  displacement  Wi  for  the  FSDPT  zone 
is  defined  in  Eq.  4.21,  and  the  displacement  W2  for  the  CPT  zone  is  defined  in  Eq.  4.12. 
These  relationships  may  be  used  along  with  Eq.  4.13  to  express  the  contribution  to  the 
potential  energy  U  of  the  spring  as 


=  2^3 


{^5,>  ^^W2> 


T  T  T  T 


(4.40) 
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This  leads  to  stiffness  subraatrix  contributions  given  by 

['^5,5,1,  =*5(05,}  K,} 


(4.41) 


=  -*5{%}  {45,} 


=  *5{%)  f%i 


where  {  a.  }  and  {  a  }  are  evaluated  at  (xi,yi)  and  (X2,y2)  respectively.  These  subraatri- 

^1  ^2 

ces  are  merged  into  the  global  stiffness  matrix  having  the  form 


1 

^1  2 

^1 3 

^1  4 

5  1 

^2  1 

Zjl  j 

^2  2 

^2  3 

^2.4, 

^2  5  1 

Zjjj  j 

^3  1 

^3  2 

JjZ] 

^33 

^34 

^3  5  1 

3j  J,  j 

^4  1 

^4  2 

^jZj 

^43 

^4  4 

^4  5  1 

^5  1 

Jjl, 

^5  2 

JjZj 

^53 

^5.4. 

^5  5  1 

^  1 
W,  1, 

K  . 

K  . 

K  . 

1 

(4.42) 


for  each  vertical  spring  connection  that  is  used 

Since  a  CPT  zone  with  a  symmetric  cross  section  (as  used  here)  does  not  have  any  in¬ 
plane  displacement  at  its  mid-surface,  there  is  no  necessity  to  enforce  in-plane  transla¬ 
tional  compatibility  with  the  FSDPT  zone.  This  leaves  only  rotational  or  slope  compati¬ 
bility  to  be  enforced.  The  rotation  61  about  the  attach  line  for  the  FSDPT  zone  is  defined 
in  Eq.  4.35,  and  the  rotation  62  about  the  attach  line  for  the  CPT  zone  is  defined  by  Eq. 
4.10.  Let  a  rotational  spring  with  stiffness  act  in  resistance  to  rotation  about  the  attach 
line  and  connect  point  (x,,yi)  on  zone  1  with  point  (x2,y2)  on  zone  2.  Substitution  of  Eqs. 
4.35  and  4.10  into  Eq.  4.37  results  in  the  potential  energy  contribution  of  the  spring  given 
in  Eq.  4.38  where  {^3  }  and  {q^  }  are  replaced  with  {q^}  ,  replaced  with 

-{a  }  ,and  {a.}  is  replaced  with  -  {a  }  .  The  resulting  stiffness  submatrix  con- 

tributions  are  given  by 
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^3  3  ^3  4 

*^4,4, 


[Jf“|  i 

sin^aj  {a3  }  {a^^  }  -costtjSinaj  {Oj  }  {a^  } 

r  nr 

-costtjsinaj  {a^}  {a-^}  cos^aj  {(^4} 


K 

K. 


3,W2 


=  kn 


T  T-[ 

sinajsina2{a3  }  -costtjsinaj  {03  } 


-cosajSina2{a4  )  J  +  cosajCosa2  {^4^}  {%^,y} 


T-[ 


(4.43) 


T  T 


]  =t, 


|^sinajsina2{a^^^}  {a^  }  -cosa2sinaj  {<23^}^ 

j^-cosajSina2{a^^j^}  {a^  }  +  cosajCosa2 {«4|}^ 


[^w,w]  =  J  {a^^y-cosa^sma^{a^^  J 

T 

-cosa^sina^ia^^J  {a^^J  +  cos2a2 


where  {a,}  and  {a.}  are  evaluated  at  (Xi,>'i),  and  {a  and  {a^  }  are  evalu- 

ated  at  (X2,y2)-  These  submatrices  are  merged  into  the  global  stiffness  matrix  having  the 
form  shown  in  Eq.  4.42  for  each  rotational  spring  connection  that  is  used.  Note  that  if  a 
rotational  spring  of  stiffness  is  used  to  resist  rotation  about  an  axis  perpendicular  to  the 
attach  line,  Eq.  4.43  may  be  used  to  find  the  stiffness  contributions  with  sina  replaced  by  - 
cosa,  cosa  replaced  by  sina,  and  k^  replaced  by  ^4. 
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4.6  Boundary  Conditions 

Boundary  conditions  along  the  wing  root  may  be  imposed  exactly  through  specific 
selection  of  displacement  polynomials,  or  approximately  through  the  use  of  computational 
springs.  For  a  CPT  zone,  exact  zero  displacement  along  the  root  line  y=0  may  be  achieved 
by  using  displacement  polynomials  that  exclude  all  terms  containing  A  cantilevered 
condition  along  the  same  line  may  be  achieved  exactly  by  using  polynomials  that  exclude 
fill  terms  containing  y°  and  y*.  For  a  FSDPT  zone,  a  cantilevered  condition  may  be 
achieved  by  using  displacement  polynomials  that  exclude  all  terms  containing  y°. 

In  using  springs,  a  zone  is  constrained  at  points  along  the  root  line  which  becomes 
the  attach  line.  At  a  particular  attach  point  on  a  FSDPT  zone,  five  different  springs  may  be 
used  to  constrain  vertical,  in-plane  translational,  and  rotational  displacements.  Assuming 
the  FSDPT  zone  to  be  zone  #1,  the  stiffness  contribution  of  a  linear  spring  with  stiffness 
resisting  vertical  displacement  is  given  by  the  first  submatrix  expression  in  Eq.  4.24.  The 
stiffness  contribution  of  a  linear  spring  with  stiffness  resisting  in-plane  translation  per¬ 
pendicular  to  the  attach  line  is  given  by  the  first  submatrix  expression  in  Eq.  4.33.  Like¬ 
wise,  the  contribution  of  a  rotational  spring  with  stiffness  k^  resisting  angular 
displacement  about  the  attach  line  is  given  by  the  first  submatrix  expression  is  Eq.  4.39. 
As  mentioned  in  Section  4.4,  springs  with  stiffness  coefficients  k^  and  k^  may  also  be  used 
respectively  to  constrain  in-plane  translations  parallel  to  the  attach  line  and  angular  dis¬ 
placements  about  an  axis  perpendicular  to  the  attach  linein  a  similar  manner.  All  of  these 
stiffness  contributions  are  merged  into  the  global  stiffness  matrix.  Note  that  a  very  large 
spring  stiffness  coefficient  will  force  displacement  to  approximately  zero.  However,  ill- 
conditioning  of  the  stiffness  matrix  results  from  using  coefficients  that  are  too  large. 


Chapter  5: 

Spar  Web  Stiffness  Approximation 


5.1  Overview 

This  chapter  discusses  the  use  of  an  equivalent  sandwich  core  to  approximate  the 
stiffness  contributions  of  an  array  of  spar  webs.  Discussed  first  is  the  case  where  spar 
webs  provide  only  transverse  shear  stiffness  to  the  wing.  Presented  next  is  the  derivation 
of  the  stiffness  contribution  of  a  sandwich  core  structure.  Finally  discussed  are  the  corre¬ 
lation  factors  allowing  the  stiffness  matrix  contribution  of  an  array  of  webs  to  be  approxi¬ 
mated  by  a  sandwich  core  over  an  associated  panel  area. 


5.2  Spar  Web  Shear  Stiffness 

Recall  from  Chapter  3  that  for  a  FSDPT  zone  the  potential  energy  contribution  due  to 
an  infinitesimal  element  of  the  jlth  spar  web  layer  is  given  by 


dUj^  =  2^y,(Tl,z) 


I  [2],/ 1  jdi]dz 

%1Z  ^  ^Tlz 


(5.1) 


The  composition  of  the  web  layer  constitutive  matrix  [U]  ji  is  given  by  (Elq.  B.166) 


2ii  2i6 
2i6  266 


(5.2) 


Now  it  is  reasonable  to  allow  the  non-shear  terms  of  [01  ji  to  be  approximately  zero  since 
spar  webs  may  contribute  very  little  to  in-plane  stiffness.  Using  this  assumption,  the  web 
is  turned  into  a  shear  web,  and  the  only  nonzero  constitutive  term  is  S66-  Thus,  the  poten¬ 
tial  energy  expression  in  Eq.  5. 1  may  be  reduced  to 


dUji  =  \tj^m,z)U66ylz^^dz 


(5.3) 


61 


Using  Eq.  2.20  for  dr\,  using  Eq.  3.51  for  and  letting  the  layer  thickness  be  a  linear 
function  of  y  gives 

(5.4) 

*V7  *VZ  \  ‘  / 


where 


[TR,,]  = 


(s^A)  (s'AcA) 
(s'AcA)  (c^A) 


Using  Eqs.  3.7, 3-9,  and  3.11-3.13  let  us  define 


[WJ  {q'} 


where 


T  T 

%  0 

T  T 

0  fl4  a^y 


T  rp  rj,  rp 

{q'}  =  {93, 


The  matrix  [W,],  2  x  iNq-i+Nq^i+Nq^)  in  dimension,  is  partitioned  into  subvectors  of  dimen¬ 
sion  1  X  Nq^  (n=3,4,5).  Substituting  Eq.  5.6  into  Eq.  5.4  gives 


dUj,  =  ^0<WS66{9'}Ow,]  ITR^]  IW,]  {q')dydz 


Integrating  over  the  area  of  the  web,  first  with  respect  to  z  and  then  with  respect  to  y,  gives 


2  cos  A 


[WJ  W)dzdy 


(5.10) 


where  the  limits  of  z  integration  are  the  depth  distributions  of  the  lower  and  upper  spar 
caps  (see  Eqs.  B.133  and  B.144).  The  z  integration  may  be  performed  analytically  (Eq. 
B.  137)  leading  to  stiffness  matrix  contribution  of  the  jlth  spar  web  layer  in  terms  of  a  spar 
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line  integral  given  by 

~h,(x.y)],.,(y)  [W/ITR,,]  lW^]dy 

(5. 11) 

where  x  along  the  spar  line  is  a  function  of  y  according  to  Eq.  2.21.  This  leads  to  linear 
combinations  of  spar  line  integrals  defined  in  Eq.  2.24.  If  a  wing  contains  many  spars,  the 
evaluation  of  web  contributions  to  the  stiffness  matrix  can  be  quite  time  consuming  when 
done  web  by  web  using  Eq.  5.11. 

5.3  Sandwich  Core  Shear  Stiffness 

Let  us  examine  the  case  where  a  sandwich  core  is  used  to  model  the  internal  structure 
of  the  wing  between  the  composite  skins  as  shown  in  Figure  5.1.  The  motivation  here  is  to 
replace  an  array  of  many  spars  and  ribs  by  a  single  equivalent  core,  whose  stiffness  can  be 
evaluated  once.  Typical  honeycomb  core  structures  theoretically  provide  only  transverse 
shear  stiffness  in  bending.  The  elastic  moduli  and  in-plane  shear  stiffness  of  the  core  are 
assumed  to  be  zero  (AL69).  Only  the  shear  moduli  and  Gy^  are  non-zero.  The  core 
constitutive  relationship  may  be  expressed  by 


Figure  5. 1:  Composite  Sandwich  Structure 
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(5.12) 


When  the  core  is  rotated  by  an  angle  (3  from  the  x-y  reference  axes  as  shown  in  Figure  5. 1 , 
this  relationship  may  be  restated  by 


where 


(  <’.■  1  = 

Q44  ^45 

|Tf„) 

045  055 

l  Y  i 

iyz 

(5.13) 


044  =  G«sin2P  +  Gy,cos2(3 

(5.14) 

045  =  (G„-Cy,)sinpcosp 

(5.15) 

Qss  =  G„cos2p  +  Gy^sin2p 

(5.16) 

(JH92).  If  the  core  is  assumed  to  be  isotropic  such  that  G„ 
becomes 


Gy,  =  G,  then  Eq.  5.13 


(5.17) 


The  constitutive  relationship  for  a  core  with  shear  stiffness  along  the  direction  of  spar 
lines,  T| ,  at  an  angle  A  from  the  y  axis  may  be  given  by 


o 


r\z 


(5.18) 


The  contribution  to  the  potential  energy  U  of  such  a  rotated  core  is  given  by 


(5.19) 


(AL69,JH92).  Using  Eqs.  3.51  and  5.5  this  may  be  rewritten  as 
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Z  V  .X 


iyz 


ITRjjl 


lyz 


dxilydz 


(5.20) 


The  z  integration  may  be  replaced  by  the  core  depth  cix,y)  shown  in  Figure  5.1.  Substitut¬ 
ing  this  and  Eq.  5.6  into  Eq.  5.20  gives 


(5.21) 


Finally,  the  resulting  stiffness  matrix  contribution  in  terms  of  an  area  integral  is  given  by 


[K]^  =  Gllcix,y)  [w/[TR^2^  [W^]dxdy 


(5.22)  . 


5.4  Spar  Web  Array  -  Sandwich  Core  Stiffness  Correlation 

Now  let  us  consider  an  array  of  spar  webs  where  all  Nwb  webs  are  equally  spaced  and 
oriented  at  the  same  angle  A  to  the  y  axis.  Referring  to  Eqs.  5.11  and  5.22,  it  seems  rea¬ 
sonable  to  assume  there  must  be  a  way  to  correlate  the  stiffness  matrix  contributions  of 
such  an  array  with  the  stiffness  matrix  contributions  of  a  sandwich  core  structure  in  an 
approximate  manner.  Doing  so  would  effectively  reduce  the  computationally  intensive 
stiffness  calculations  of  several  spar  webs  composed  of  multiple  layers  to  relatively  sim¬ 
ple  calculations  over  a  small  number  of  single-layered  panels. 

Using  Eq.  5.11,  an  averaged  spar  web  stiffness  matrix  can  be  determined  for  an  array 
of  spar  webs  corresponding  to  a  particular  trapezoidal  panel  area.  The  web  “density  fac¬ 
tor”  d  may  be  determined  from  the  number  of  webs  per  chord  dimension  of  the  panel. 
Referring  to  Figure  1.3,  this  is  given  by 

d  = 

X^L  ~  ^FL  (5.23) 

We  assume  that  all  spar  webs  have  NU  layers  each  with  the  same  thickness.  Therefore,  an 
average  web  thickness  t„b(y)  can  be  determined  by  summing  over  the  layer  thicknesses  of 
one  web.  We  allow  for  each  composite  layer  to  have  a  different  shear  stiffness  due  to  prin¬ 
cipal  material  axis  orientation.  An  average  shear  stiffness  for  the  whole  web  can  be  deter¬ 
mined  at  any  y  coordinate  along  one  web  using 
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‘wfcV// — J  '  '  (5.24) 

Using  these  averaged  values,  the  averaged  spar  web  stiffness  matrix  is  given  by 

d  y,  \  r  t  /  / 


Thus,  Eq.  5.22  becomes  approximately  equivalent  to  Eq.  5.25  when  G  =  Q^y  and 


(5.25) 


(5.26) 


Therefore,  the  stiffness  matrix  contribution  of  an  approximately  equivalent  core  over  the 
panel  area  associated  with  the  web  array  is  given  by 


[hu(x.y)  -/iiU.j’)]  {V/,]dxdy 


yx 


(5.21) 


A  whole  array  of  evenly  distributed  spar  webs  having  the  same  web  construction  can 
now  be  replaced  by  a  single  layer  of  equivalent  core  material.  Evaluation  of  the  core  stiff¬ 
ness  contribution  is  done  using  the  same  integral  tables  used  to  evaluate  contributions  of 
cover  skins,  thus  leading  to  substantial  saving  of  computational  resources.  A  similar  sim¬ 
plification  may  be  made  for  an  array  of  rib  webs. 


Chapter  6: 
Test  Case  Results 


6.1  Overview 

This  chapter  presents  numerical  results  obtained  with  the  new  equivalent  plate  model¬ 
ing  capabilities  focusing  on  two  different  test  cases.  The  first  test  case  involves  a  swept, 
thick,  high  aspect  ratio  wing  of  the  type  used  for  subsonic  transports.  This  case  was  used 
to  study  the  effects  of  using  Ritz  polynomials  of  different  orders  for  different  displacement 
fields  and  to  study  the  accuracy  of  replacing  detailed  spar  web  modeling  by  an  equivalent 
core  are  studied.  A  Boeing  HSCT  candidate  wing  represents  a  configuration  made  of  a 
low  aspect  ratio  inner  section  and  a  high  aspect  ratio  outer  wing,  densely  packed  with  sup¬ 
porting  internal  spars.  Effectiveness  and  accuracy  of  equivalent  core  modeling  are  tested 
with  this  wing. 


Figure  6. 1:  Turner-Martin-Weikel  Wing 
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6.2  Swept,  Thick,  High  Aspect  Ratio  Wing 

A  all  aluminum  subsonic  transport  type  wing,  for  which  experimental  and  numerical 
results  are  available,  is  discussed  in  TMW64  and  LIV194  .  Shown  in  Figure  6.1  is  the 
geometry  of  this  wing  which  has  a  symmetric  airfoil  (see  Appendix  D  for  more  detail). 
This  wing  is  used  here  first  to  test  the  effect  of  varying  the  order  of  polynomials  used  for 
different  displacement  and  rotation  fields  in  the  FSDPT  model.  Recall  from  Chapter  3  that 
there  are  3  independent  linear  mid-plane  displacements  Uq,  Vq,  Wq  and  2  independent  rota¬ 
tional  displacements  Xj/;,,  X|/j,  in  the  model.  Obviously,  the  vertical  deflection,  stress,  and 
natural  frequency  predictions  of  the  model  will  deteriorate  if  polynomials  of  low  order  are 
used.  However,  the  question  is  whether  it  is  possible  to  use  lower  order  polynomials  for 
the  shear  rotation  fields  or  any  other  deformation  field  without  severely  degrading  the 
overall  accuracy.  If  this  is  so,  the  overall  order  of  wing  models  may  be  reduced  leading  to 
computational  savings  in  storage  and  CPU  time. 

The  numerical  tests  involved  first  varying  the  powers  of  the  polynomials  used  for  \|/^ 
and  Yj,  alone.  Then  the  powers  of  the  polynomials  used  for  Mq  and  Vq  alone  were  varied. 


Figure  6.2:  Vertical  Displacement  along  Rear  Spar  with  Varying 
Orders  and\\fy  Rotation  Field  Polynomials 
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The  reference  case  was  based  on  using  complete  6th  order  polynomials  for  all  5  displace¬ 
ment  fields.  Root  cantilever  boundary  conditions  along  y=0  were  enforced  by  omitting 
terms  containing  y®  from  the  displacement  polynomials,  thus  leaving  21  terms  for  each 
displacement  for  a  total  of  105  degrees  of  freedom.  A  1  lb.  vertical  load  was  applied  at  the 
outboard  point  on  the  trailing  edge. 

Available  experimental  test  data  and  finite  element  results  are  used  here  for  compari¬ 
son.  Figure  6.2  shows  the  vertical  deflections  along  the  rear  spar  with  rotation  \\f^  and 
polynomials  varying  in  order  from  second  to  sixth  order.  Shown  in  Figure  6.3  is  the  nor¬ 
mal  stress  in  the  primary  spar  direction  along  the  root  chord  from  front  to  rear. 

Note  that  the  finite  element  stress  results  do  not  correlate  well  in  the  region  of  higher  stress 
gradient  due  to  choice  of  mesh  size.  Natural  frequency  results  are  given  in  Table  6.1.  The 
results  show  that  although  deflection  and  natural  frequency  correlation  is  still  good  for 
polynomials  of  4th  order  (10  terms)  stress  correlation  deteriorates  quickly  in  the  area  of 
high  stress  gradients  at  the  inboard  trailing  edge.  Still,  with  4th  order  rotation  field  poly¬ 
nomials,  overall  stress  accuracy  is  comparable  to  the  finite  element  results. 
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x/(root  chord) 


Figure  6.3:  Skin  Stress  along  Root  Chord  with  Varying 
Orders  and\yy  Rotation  Field  Polynomials 
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Table  6.1:  Natural  Frequencies  (in  Hz)  with  Varying  Orders  of 
and  Y),  Rotation  Field  Polynomials 


Frequency  # 
(Shape) 

FEM 

6th 

Order 

5th 

Order 

4th 

Order 

3rd 

Order 

2nd 

Order 

#1  (1st  Bending) 

115.6 

114.7 

115.1 

115.8 

118.1 

125.7 

#2  (Fore/  aft) 

317.6 

312.4 

312.4 

312.4 

312.4 

312.4 

#3  (1st  Torsion) 

418.4 

428.9 

430.3 

431.9 

441.6 

484.8 

#4  (2nd  Bending) 

576.4 

575.3 

576.9 

582.8 

591.9 

715.0 

#5  (2nd  Torsion) 

1086 

1125 

1136 

1157 

1251 

1251 

When  the  order  of  polynomials  used  for  Uq  and  Vq  was  changed,  no  change  from  the 
reference  case  was  observed  in  vertical  deflections  or  stresses  using  polynomials  down  to 
1st  order  (1  term).  The  2nd  natural  frequency,  corresponding  to  fore/aft  motion,  was  the 
only  value  that  changed  at  all  throughout  the  comparison.  This  data  is  shown  in  Table  6.2. 
Such  results  could  be  expected  because  of  the  symmetry  of  airfoil  cross  section.  If  inaccu¬ 
rate  prediction  of  the  fore/aft  vibration  frequeney  can  be  tolerated  (whieh  is  indeed  the 
case  when  classical  flutter  is  involved),  then  only  53  degrees  of  freedom  (21  for  Wq,  15  for 
Y;t  and  y>»  1  for  and  Vo)  are  needed  for  the  wing.  This  is  approximately  a  50% 

reduction  in  model  size  (compared  with  the  5  x  21  =  105  dof  model)  leading  to  just  a  small 
reduction  in  accuracy. 


Table  6.2:  Second  Natural  Frequency  (in  Hz)  for  Varying 
Orders  of  Mq  and  Vo  Displacement  Polynomials 


6th 

Order 

5th 

Order 

4th 

Order 

3rd 

Order 

2nd 

Order 

1st 

Order 

312.4 

316.3 

329.7 

364.1 

428.9 

428.9 
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Figure  6.4:  Vertical  Displacement  along  Rear  Spar 
Using  an  Equivalent  Shear  Core 

The  Tumer-Martin-Weikel  wing  is  also  used  here  as  a  test  case  to  assess  the  accuracy 
of  an  equivalent  core  representation  for  shear  web  stiffness.  Stiffness  of  the  array  of  5 
spar  webs  parallel  to  the  leading  edge  is  approximated  by  a  single  shear  core  using  Eq. 
5.27.  The  same  loading  condition  and  displacement  polynomials  from  the  reference  case 
above  are  used  here.  Vertical  deflections  along  the  rear  spar  are  shown  in  Figure  6.4.  Cor¬ 
relation  with  the  reference  case  is  reasonable  with  an  error  of  4.4%  at  the  outboard  point. 
Stress  results  along  the  root  chord  are  shown  in  Figure  6.5.  Correlation  is  reasonable  here 
as  well  with  an  error  of  6. 1  %  at  the  inboard  point  on  the  trailing  edge.  Natural  frequency 
results  are  reported  in  Table  6.3.  The  most  serious  discrepancies  with  respect  to  the  refer¬ 
ence  case  occur  in  the  3rd  frequency  (9.4%  error)  and  the  5th  frequency  (6.4%  error),  both 
of  which  correspond  to  torsion  modes.  However,  notice  that  these  results  have  an  error  of 
only  7.1%  and  3.0%  respectively  when  compared  to  the  finite  element  standard.  Results 
from  the  test  indicate  that  the  equivalent  core  approximation  works  well  for  this  thick, 
high  aspect  ratio  wing. 


Normal  Stress  in  Elastic 
Axis  Direction  (psi) 


15 


Figure  6.5:  Skin  Stress  aiong  Root  Chord 
Using  an  Equivalent  Shear  Core 


Table  6.3:  Natural  Frequencies  (in  Hz)  Using 
an  Equivalent  Shear  Core 


Frequency  # 
(Shape) 

FEM 

FSDPT 

Reference 

FSDPT 
Equiv.  Core 

#1  (1st  Bending) 

115.6 

114.7 

114.4 

#2  (Fore/  aft) 

317.6 

312.4 

312.4 

#3  (1st  Torsion) 

418.4 

428.9 

388.7 

#4  (2nd  Bending) 

576.4 

575.3 

568.9 

#5  (2nd  Torsion) 

1086 

1125 

1053 
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6.3  HSCTWing 

A  simplified  model  of  a  representative  HSCT  wing  was  used  to  assess  the  accuracy  of 
CPT  based  modeling  of  wings  (LSB93).  The  simplified  wing  is  symmetric  about  the  mid¬ 
plane  and  has  a  linearly  varying  depth  from  root  to  tip.  Its  planform  layout  is  shown  in 
Figure  6.6  (see  Appendix  D).  Previous  HSCT  numerical  tests  included  finite  element 
results  obtained  using  the  airframe  structural  optimization  code  ELFINI  (LSB93,LP87). 
With  42  spar  webs,  this  wing  provides  a  challenging  test  case  for  the  equivalent  core 
approximation.  If  the  preparation  of  stiffness  contributions  due  to  42  webs,  done  one  by 
one  as  in  LSB93,  can  be  replaced  by  the  evaluation  of  stiffne.ss  contribution  of  one  core, 
then  CPU  time  savings  can  be  significant. 

Figure  6.7  shows  vertical  displacements  for  the  HSCT  wing  along  lines  of  constant  y. 
The  results  show  that  u.se  of  the  equivalent  core  gives  good  correlation  with  displacements 
based  on  individual  spar  webs.  Using  the  equivalent  core  leads  to  a  discrepancy  of  only 
2.9%  at  the  outer  wing  tip.  Stress  results  are  given  in  Figure  6.8  for  the  skin  layer  ori¬ 
ented  at  0°  to  the  inboard  wing  spar  direction.  Use  of  the  equivalent  core  gives  good  cor¬ 
relation  here  as  well. 


1200  1400  1600  1800  2000  2200  2400  2600  2800 

X  (in)  -  chord  direction 


Figure  6.6:  HSCT  Wing  Planform 
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Figure  6.7:  HSCT  Wing  Vertical  Displacements 


Table  6.4:  HSCT  Natural  Frequencies  (in  Hz) 


Frequency  # 

ELFINI 

FSDPT 

Reference 

FSDPT 
Equiv.  Core 

#1 

9.50 

9.50 

9.37 

#2 

23.9 

24.2 

24.0 

#3 

26.9 

24.8 

24.7 

#4 

33.7 

32.3 

32.2 

#5 

36.1 

38.8 

38.8 

#6 

42.7 

43.8 

43.0 

#7 

45.2 

50.0 

48.8 

#8 

49.6 

55.5 

55.5 
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Figure  6.8:  Normal  Stress  in  Layer  Oriented  at  0° 
to  Inboard  Wing  Spar  Direction 

Table  6.4  shows  a  comparison  of  natural  frequencies  for  ELFINI,  FSDPT  with  normal 
spar  webs,  and  FSDPT  with  the  equivalent  shear  core.  The  frequencies  are  dependent 
upon  the  stiffness  of  the  springs  used  to  constrain  displacements  at  the  wing  root.  Using 
the  spring  stiffness  coefficients  listed  in  Appendix  D,  the  FSDPT  capability  u.sing  spar 
webs  matches  the  ELFINI  capability  reasonably.  The  equivalent  core  results  correlate 


Table  6.5:  Stiffness  Matrix  Assembly  Time  (CPU  seconds) 


#  of  terms  for 
Displacements 

42  Spar  Webs 

Equivalent  Core 

10 

2893.5 

14.1 

15 

1861.5 

10.8 

21 

765.7 

3.3 

75 


well  with  both  of  these.  Both  tests  with  the  Tumer-Martin-Weikel  wing  and  HSCT  wing 
demonstrate  the  accuracy  of  modeling  when  an  array  of  spar  webs  is  replaced  by  an  equiv¬ 
alent  core. 

The  effectiveness  of  using  an  equivalent  core  is  shown  in  a  comparison  of  the  com¬ 
puter  resources  required  to  assemble  the  stiffness  matrix  contributions  for  an  array  of  spar 
webs  and  an  equivalent  core.  Computations  were  performed  on  an  HP  Apollo  700  work¬ 
station.  Tables  6.5  compares  the  CPU  time  used  in  assembling  the  stiffness  matrix  for  42 
spar  webs  vs.  the  time  used  in  assembling  the  stiffness  matrix  for  an  equivalent  core.  As 
hoped  for,  the  computational  savings  are  significant. 


Chapter  7: 
Conclusion 


7.1  Summary 

This  work  has  presented  recent  improvements  in  equivalent  plate  modeling  of  aircraft 
wings  for  use  in  multi-disciplinary  design  optimization.  A  mathematical  wing  model  was 
outlined  using  energy  principles  and  the  Ritz  solution  method.  Formulations  of  mass, 
stiffness,  and  load  contributions  were  given  using  both  Classical  Plate  Theory  and  First 
Order  Shear  Deformation  Plate  Theory.  These  formulations  allow  each  displacement  field 
in  the  model  to  be  independently  approximated  by  a  polynomial  Ritz  series  of  a  chosen 
order.  Additionally,  a  formulation  for  determining  the  spring  stiffness  contributions  from 
the  connection  of  a  CPT  zone  to  a  FSDPT  zone  was  derived.  A  method  for  approximating 
the  stiffness  of  an  array  of  spar  webs  using  an  equivalent  core  was  also  developed. 

Results  are  presented  from  numerical  tests  performed  on  two  different  wings.  Tests 
of  the  Turner-Martin-Weikel  wing  showed  that  low  order  polynomials  may  be  used  for 
linear  in-plane  displacements  at  the  mid-surface,  resulting  in  retained  accuracy  from  a 
model  of  significantly  reduced  order.  This  wing  was  also  used  to  show  that  an  equivalent 
core  can  be  used  in  place  of  an  array  of  spar  webs  with  acceptable  accuracy.  Tests  of  a 
simplified  HSCT  wing  validated  the  accuracy  and  demonstrated  the  computational  effi¬ 
ciency  of  the  equivalent  core  approximation  applied  to  a  wing  with  many  spar  webs. 

The  findings  here  serve  to  advance  the  capabilities  of  equivalent  plate  wing  modeling. 
These  new  developments  allow  wings  of  increasing  complexity  to  be  modeled  with  a  bet¬ 
ter  combination  of  accuracy  and  efficiency  than  could  be  achieved  previously.  Ultimately, 
this  builds  support  for  the  use  of  equivalent  plate  modeling  as  the  primary  wing  structural 
analysis  tool  in  preliminary  design  optimization. 
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7.2  Future  Work 

Future  extension  of  this  work  should  include  thorough  numerical  testing  of  the  formu¬ 
lation  for  spring  connection  of  FSDPT  and  CPT  zones.  This  new  capability  should  be 
tested  on  several  wing  configurations  with  available  data  for  comparison.  Also,  for  this 
capability  to  be  used  in  design  optimization,  analytic  sensitivities  with  respect  to  the 
design  variables  must  be  obtained  and  implemented  for  any  new  parameters  introduced 
into  the  model.  The  developed  code  must  also  be  modified  to  allow  it  to  be  linked  with 
optimizing  routines. 
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Appendix  A: 

Analytical  Foundations  of  an  Equivalent 
Plate  Model  Using  CPT 


A.1  Overview 

In  this  appendix  are  derived  the  equations  necessary  to  construct  mass  and  stiffness 
matrices  for  wing  structural  components  based  on  Classical  Plate  Theory.  This  appendix 
starts  with  a  discussion  of  the  assumptions  of  Classical  Plate  Theory  and  the  ensuing 
strain-displacement  relationships.  Following  this  foundation  are  the  derivations  of  the 
mass  and  stiffness  matrix  equations  for  a  single  skin  layer.  Next  are  the  derivations  of  the 
mass  and  stiffness  matrix  equations  for  a  single  spar.  Presented  last  are  the  derivations  of 
the  mass  and  stiffness  matrix  equations  for  a  single  rib. 


A.2  Assumptions 

Classical  Plate  Theory  (CPT)  is  the  name  given  to  the  small-deflection  theory  of 
bending  of  elastic  thin  plates  (RE84).  The  first  assumption  of  CPT  for  a  plate  in  the  x-y 
reference  plane  is  that  the  displacements  of  the  mid-surface  are  small  compared  to  the 
thickness  of  the  plate,  thus  meaning  the  slope  of  the  deflected  surface  is  very  small.  The 
second  assumption  is  that  the  plate  mid-surface  remains  unstrained  after  bending.  The 
third  assumption,  known  as  the  Kirchoff-Love  assumption,  is  that  planes  normal  to  the 
plate’s  mid-surface  will  remain  plane  and  normal  to  the  plate’s  mid-surface  after  the  plate 
is  deformed  (RE84,J075).  This  assumption  is  equivalent  to  ignoring  the  shear  strains  y„ 
and  Yy^.  The  fourth  assumption  is  that  the  transverse  normal  strain  is  negligible  since 
the  transverse  loads  are  transferred  primarily  to  bending  strains.  In  other  words,  the  trans¬ 
verse  deflection  does  not  vary  across  the  plate  thickness.  The  final  assumption  is  that  the 
transverse  normal  stress  (J^  may  be  neglected  since  it  is  small  compared  to  the  other 
stresses. 

Certain  assumptions  must  be  made  about  the  wing  being  modeled  as  well.  The  first  of 
these  assumptions  is  that  no  spar  webs  or  rib  webs  will  be  included.  Only  skins,  spar  caps, 
and  rib  caps,  will  be  accounted  for.  See  Figure  1.1.  The  second  assumption  is  that  the 
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wing  is  thin  and  symmetric  with  respect  to  the  x-y  plane.  The  final  assumption  is  that  the 
wing  skins  are  thin  with  respect  to  the  wing  depth  allowing  all  skin  layers  in  each  of  the 
upper  and  lower  skins  to  have  the  same  depth  distribution  (LIV90). 


A.3  Strain-Displacement  Relationships 

The  established  assumptions  allow  the  plate  problem  to  be  considered  under  plane 
stress  conditions  where  the  in-plane  strains  £yy>  and  are  the  only  strains  of  concern. 
It  may  also  be  concluded  that  the  normal  displacement  in  the  z  direction  is  a  function  of 
only  X  and  y;  therefore 


w  -  w  {x,  y) 


(A.1) 


Because  of  the  Kirchoff-Love  assumption  and  the  symmetry  of  the  wing,  the  inplane  dis¬ 
placements  u  and  V  are  given  as 


u 


-zw 


'X 


dw 


(A.2) 

(A.3) 


The  strains  are  then  given  as 


e 


XX 


du 

dx 


=  -zw. 


XX 


(A.4) 


e 


yy 


dv 

dy 


=  -zw,. 


yy 


(A.5) 


_  du  dv 
^■*3'  dx 


-2zw, 


xy 


(A.6) 


It  is  evident  that  all  the  displacements  and  strains  are  functions  of  the  independent  trans¬ 
verse  displacement  w(x,y).  These  relationships  form  the  basis  for  the  stiffness  and  mass 
matrix  formulations  (LIV90). 
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A.4  Displacement  Function 

To  use  the  Ritz  method  described  in  Chapter  1,  it  is  necessary  to  define  the  displace¬ 
ment  function  w(x,y,t)  as  a  polynomial  series  having  terms  using  Eq.  1.9.  Thus, 


—  mQi  nq- 

wix,y,t)  =  ^qy^(t)i‘X  y 

i=l 


(A.7) 


where  the  time  dependent  coefficients  ^„,(t);  become  the  generalized  displacements  solved 
for  in  the  Ritz  formulation,  and  the  powers  m^;  and  are  predetermined  (LIV90).  This 
may  be  expressed  more  simply  in  vector  form  as 


T 

wix,y,t)  =  {a^(x,y)}  {q^it)} 


(A.8) 


where 

T  mq,  nq, 

=  {...,x  y  ,...} 

(A.9) 

and  {qjf)}  is  a  column  vector.  Both  vectors  contain  Nq^  terms.  This  vector  form  for 
w{x,y,t)  will  be  used  in  the  mass  and  stiffness  matrix  formulations. 


A.5  Mass  Matrix  Contributions  of  a  Skin  Layer 


Neglecting  rotary  inertia  the  contribution  to  the  kinetic  energy  7  of  an  infinitesimal 
skin  element  dx  by  dy  of  the  y/th  skin  layer  is 


dTji  =  :^pjitji{x,y)w  {x,y,t)dxdy 


(A.10) 


where  is  the  constant  material  density  of  the  skin  layer,  and  tji{x,y)  is  the  thickness  of  the 
skin  layer  given  in  Eq.  1.2  as 


m., 


r,,0.y)  =  ‘y 


^  =  1 


(A.11) 


(LIV90).  Substituting  Eq.  A.8  into  Eq.  A.  10  gives 
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1  ^ 

dTji  =  -Pjitj^{x,y)  {qj  {aj  {aj  {qjdxdy 
Integrating  this  over  the  whole  trapezoidal  panel  area  yields 

T  T 

Tji  =  \ 9ji\ I tji  (X,  y)  { qj  {aj{aj  { qj  dxdy 

yx 


(A.12) 


(A.13) 


Equivalently, 


yx 


(A.14) 


where  co  is  the  natural  vibrational  frequency  of  the  skin  layer.  The  variational  extremum 
condition  from  Eq.  1.10  applied  here  becomes 


=  0 


(A.15) 


yielding  the  mass  matrix  for  the  jlth  skin  layer  given  as 


^^Skhi  =  9ji\\tji^x^y)  {aj  dxdy 


yx 


(A.16) 


Substitution  of  Eq.  A.  1 1  for  tji(x,y)  gives 


k=i  yx 


(A.17) 


The  mass  matrix  term  in  the  ith  row  and  jth  column  for  the  jlth  skin  layer  may  be 
expressed  as 


where 


k=i  yx 


(A.18) 


m  =  mq.  +  mq-  +  mr^ 


(A.19) 
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n  =  nq-  +  nqj  +  nti^ 

Let  us  define  a  family  of  integrals  over  the  area  of  the  trapezoidal  panel  where 

=  ^^x”'y'dxdy 

yx 

Then  the  final  expression  for  the  ijth  term  of  is 

Ntj, 


(A.20) 


(A.21) 


k=l 


(A.22) 


Appendix  C  discusses  how  the  family  of  integrals  Ij-[^{m,n)  is  analytically  derived.  The 
mass  matrix  for  all  the  skin  layers  combined  is  obtained  by  summing  the  contribu¬ 

tions  of  all  the  layers  given  by  Eq.  A.22.  Because  the  wing  is  assumed  to  be  symmetric 
with  respect  to  the  jc-y  plane,  the  total  skin  contribution  can  be  found  by  summing  the  con¬ 
tributions  from  the  layers  of  the  upper  skin  and  then  multiplying  by  2.  This  matrix  is  then 
appropriately  merged  into  the  global  mass  matrix 


A.6  Stiffness  Matrix  Contributions  of  a  Skin  Layer 

The  contribution  to  the  potential  or  strain  energy  t/  of  an  infinitesimal  skin  element  dx 
by  dy  of  the  jlth  skin  layer  is 


1 


dUji  =  ^{K}  [D]ji{K}dxdy 


(A.23) 


where 


w, 


XX 


{K}  =  -1 

2w, 


(A.24) 


’.cy 


using  Eqs.  A.4,  A.5,  and  A.6,  and  [D]jiis  the  plate  bending  stiffness  matrix  defined  by 

[D] =  I lQ]/dz 


(A.25) 
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(J075,  LIV90).  The  jlth  layer’s  constitutive  matrix  [Q]p  as  discussed  in  Section  1.5,  has 
the  form  r  n 


Qn  Q\i  ^16 
Q\i  Qll  ^26 


2i6  ^26  ^66 


(A.26) 


where  its  components  are  determined  for  the  x-y  axis  system  from  [Q\  ji ,  the  y/th  layer’s 
constitutive  matrix  referenced  to  its  own  principal  material  axes.  Assuming  the  jlth  layer 
to  be  an  orthotropic  lamina,  its  invariant  properties  may  be  used  to  find  [Qlyvwhen  the  prin¬ 
cipal  material  axes  are  oriented  at  some  angle  p  to  the  x-y  axes  as  shown  in  Figure  A.  1 
(J075).  The  components  of  [Q]j[  may  be  written  as 


(2ii  =  Cj  +  C2COS2P  +  C3COS4P 
Qi2  =  C4-C3COS4P 
Q22  ==  Cj  -  C2COS2P  +  C3COS4P 
Qi6  =  -^C2sin2p-C3sin4p 
Q26  =  -  ^C2sin2p  +  C3sin4p 
Qee  =  C5-C3COS4P 


(A.27) 

(A.28) 

(A.29) 

(A.30) 

(A.31) 

(A.32) 


Figure  A.  1:  Positive  Rotation  of  Principal  Material  Axes  from  x-y  Axes 
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for  which  the  invariants  are  given  as 

3Gii  +  3022  +  2012  +  4066 


C, 


8 


(A.33) 


C2  - 


Q\\  -  Q22 


(A.34) 


C3  = 


^11  +  Q22~^Qi2~^Q66 

8 


(A.35) 


^^4  = 


C5  = 


611  +  622  +  ^Qi2~  ^Qee 
8 


611  622  ~  2^12  +  ^666 
8 


(A.36) 


(A.37) 


Referring  to  Figure  A.2,  the  jlth  skin  layer  lies  between  the  coordinates  z=h/2-tji/2  and 
z=hl2+tji/2  such  that  the  z  integral  from  Eq.  A.25  is  given  by 


f(^'*''^')/2^2  ,  -  •j_f 


h 


L  3V/iv'J 


(A.38) 


Since  it  is  assumed  that  the  wing  thickness  is  much  smaller  than  the  wing  depth  (tj/h  «1), 
the  right  hand  side  of  Eq.  A.38  can  be  simplified  to  h^tj/A  (LIV90).  Substituting  Eqs.  A.24 
and  A.25,  and  the  simplified  Eq.  A.38  into  Eq.  A.23  then  gives 


dUj, 


If  hix,y)' 


T 

^xx 

’JCJC 

'  ^Q\ji 

-  2w,^^  . 

J 

dxdy 


(A.39) 


^xy 


Substituting  Eq.  A.8  results  in 
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(A.43) 
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yielding  the  stiffness  matrix  for  the  y/th  skin  layer  given  as 


=  jj 


(A.44) 


Substitution  of  Eq.  A.  1 1  for  tu{x,y)  and  the  following 


n/i,^ 


AUy)  =  X  ^  y 


(A.45) 


for  h{x,y)  (See  Eq.  1.1)  into  Eq.  A.43  yields 

Nh  Nh  Ntj, 

1^=1  jk=lk=l 

mh,^  +  mhj^  +  mt^  nhi^+nhj^  +  nt/^ 


[W]  [Q]ji[W]dxdy 


(A.46) 


Let  the  r,/th  term  of  the  matrix  [W]  (3  x  A^^in  dimension)  be 


rii//  -M  TI7/  mW(r,i)  nW(r,i) 

[W(r,  0]  =  Wir,i)x  y 


(A.47) 


Then  the  ijth  term  of  the  triple  matrix  product  [W]^[(2]j7[W]  is 


w/  •\{irr  ^  n  i  \  "JlVlr.i)  +mW(,yJ)  wW(r,0  +wW(i'.7) 

term(i,j)  =  T  \  W  (r,  i)W  {s,j)  Q.^ir,  s)  x  y 

r=i.=  i 


Upon  substitution  of  Eq.  A.48  into  Eq.  A.46,  the  stiffness  matrix  term  in  the  ith  row  and 
yth  column  for  the  j7th  skin  layer  may  be  expressed  as 


Nh  Nh  Ntj,  3  3 


ik=i  jk-lk=l  r=l 


(A.49) 


•  W  (r,  i)  W  (sj)  Qji  (r,  s)  JJx'”y"rfxJy 


where 
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m  =  mW{r,i)  +  mW  {s,j)  +  mh-i^  + mhji^  + rntf^ 


(A.50) 


n  =  nW (r,  i)  +nW{s,j)  +  nh-i^  +  nhj^.  + 


(A.51) 


Noting  that  Eq.  A.49  contains  an  integral  from  the  family  lTR{m,n)  defined  in  Eq.  A.21,  the 
final  expression  for  the  I'jth  term  of  is 


Nh  Nh  Ntj,  3  3 

ik=ljk=ik=l  r=ls=l  (A.52) 

•  W  (r,  i)  WisJ)  Qji  (r,  s)  (m,  n) 

The  stiffness  matrix  for  all  the  skin  layers  combined  is  obtained  by  summing  the 

contributions  of  all  the  layers  given  by  Eq.  A.52.  Because  the  wing  is  assumed  to  be  sym¬ 
metric  with  respect  to  the  x-y  plane,  the  total  skin  contribution  can  be  found  by  summing 
the  contributions  from  the  layers  of  the  upper  skin  and  then  multiplying  by  2.  This  matrix 
is  then  appropriately  merged  into  the  global  mass  matrix 


A.7  Mass  Matrix  Contributions  of  a  Spar 

Neglecting  rotary  inertia  the  contribution  to  the  kinetic  energy  T  of  an  element  of 
length  dr\  of  the  jsth  spar  is 

JS  2  (A  53) 

where  Tj  is  the  coordinate  along  the  spar  axis  rotated  from  the  y  axis  by  an  angle  A,  is 
the  constant  material  density  of  the  spar,  and  the  cap  area  A  is  expressed  as  a  function 
of  q  (LIV90).  Referring  to  Figure  A.3,  all  q  dependence  of  this  equation  may  be  changed 
to  y  dependence  using 

q  ^  y-^yt 

L  syj^-sy^  (A.54) 


from  which 


Figure  A. 3:  Spar  Geometry 


dr{  =dy 


L  _  dy 
cos  a 


(A.55) 


The  spar  cap  area  may  be  expressed  as  a  linear  function  of  y  by 


(A.56) 


as  shown  in  Eq.  1.3.  The  variable  x  may  also  be  expressed  as  a  function  of  y  in  the  form 

X  (y)  =  51y  +  S2 


(A.57) 


where 


51  = 


SXj^  -  sx^ 

^yR-^yi 


(A.58) 


52  = 


sXpSyR-sxRsy^ 

^yR-^yi 


(A.59) 


Combining  Eq.  A.57  with  Eq.  A.7  allows  the  displacement  w{x,y,t)  to  be  written  as 


w 


Nq. 

i  =  1 


(A.60) 


from  which  Eq.  A.9  becomes 
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Using  Eqs.  A.8,  A.55,  and  A.61,  the  expression  from  Eq.  A.53  now  becomes 

O’).,,,  {««■>  ‘  {o*}  (“».} '  t«»)  ‘‘y 

in  terms  of  Integrating  over  the  length  of  the  spar  gives 


(A.61) 


(A.62) 


(A.63) 


Equivalently, 


T,,  =  Wp,: 


,  -  sy,  Jjy=^yL 


{?w}  {aj{a\  {qjdy 


where  co  is  the  natural  vibrational  frequency  of  the  spar.  Application  here  of  the  varia¬ 
tional  extremum  condition  from  Eq.  1.10  yields  the  mass  matrix  for  they^th  spar  given  as 

I  X  r 


(A.65) 


Substituting  Eq.  A.56  gives 


[".ply.  =  Py.(;;;VJn^^»pOp-*'''ppl/l  {aj  {afdy 


(A.66) 


The  mass  matrix  term  in  the  ith  row  and  yth  column  for  the  y'^th  spar  may  be  expressed  as 


IM  (i.J)].  =  p,rA,„„f— ^ln'(Sly  +  S2)'”/.;y 

L  sp^^J^^js  ^P^js\syj^-sy^Jjy=syL 

L 


+A, 


(A.67) 


hs\  syjf  -  sy^Jh^^yi 


ri‘(Sly.S2r/''dy 


where 


m  -  mq-  +  mqj 
n  =  n^-  +  nqj 


(A.68) 

(A.69) 
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Let  us  define  a  family  of  integrals  where 

L 


I^p{m,n)  = 


syp  -  sypjy=^yL 


r;?(5iy+52)V^y 


(A.70) 


Then  the  final  expression  for  the  ijth  term  of  [M^p]p  is 

(ij)  ] =  Pis  ^\voJsp  «)  +  \piJsp  (m,  n  +  1)  ] 


(A.71) 


(LIV90).  Appendix  C  discusses  how  the  family  of  integrals  Isp(ni,n)  is  analytically 
derived.  The  mass  matrix  for  all  the  spars  combined  [M,p\,„^is  obtained  by  summing  the 
contributions  of  all  the  spars  given  by  Eq.  A.71.  Because  the  wing  is  assumed  to  be  sym¬ 
metric  with  respect  to  the  x-y  plane,  the  total  spar  contribution  can  be  found  by  summing 
the  contributions  from  the  spars  on  the  upper  wing  surface  and  then  multiplying  by  2. 
This  matrix  is  then  appropriately  merged  into  the  global  mass  matrix 


A.8  Stiffness  Matrix  Contributions  of  a  Spar 

The  contribution  to  the  potential  or  strain  energy  (/  of  an  element  of  length  dT\  of  the 
jsth  spar  is 


d'Jj.  = 


(A.72) 


where  the  longitudinal  modulus  of  elasticity,  the  cap  area  A  is  a  linear  function  of 
T],  and  h{r[)  is  the  wing  depth  distribution  along  the  spar  line  (LrV90).  It  is  necessary  to 
replace  Tj  dependence  with  y  dependence  in  this  expression.  Referring  to  Figure  A.3  and 
Eq.  A.55  the  transformation  for  the  bending  strain  is  given  by 


f syp  -  Y 

=  COS^A  -  - J - J 


(A.73) 


The  depth  distribution  may  be  written  as  a  function  of  y  by  substituting  the  spar  line  equa¬ 
tion  from  Eq.  A.57  into  Eq.  A.45  giving 


Nh 


mhik 


hiy)  =  2  iSly  +  S2)  ‘^y 


ik  =  1 


(A.74) 
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The  cap  area  is  also  easily  expressed  as  a  linear  function  of  y  as  shown  in  Eq.  A.56.  Sub¬ 
stituting  Eqs.  A.55  and  A73  into  Eq.  A.72  gives 


(A.75) 


Using  Eqs.  A.8  and  A.61  to  define  wix,y,t),  Eq.  A.75  becomes 


(A.76) 


Integrating  this  over  the  length  of  the  spar  gives 

2^ A  L  )  h=syA  ^Pjs  4 


(A.77) 


{qj  {a^yy}  {qjdy 

Application  here  of  the  variational  extremum  condition  from  Eq.  1.10  yields  the  stiffness 
matrix  for  the  ysth  spar  given  as 


IV,  -  A-^}  ‘“-U  (a.78) 


Substituting  Eq.  A.56  for  A  and  Eq.  A.74  for  h(y)  here  gives 

rjs 


ik=\  jk-l 


(A.79) 


(Sly +  52)  y  {a^,yy}  {a^,yy}  dy 


Now  differentiating  Eq.  A.61  twice  with  respect  to  y  gives 
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mqi-2  nqi 


(Sly  +  Sl)  y 


mO:  -  1  nq:  -  1 

+2mq.nq.SliSly  +  S2)  y 

mo;  nq-  -  2 

+nq.inq.-l)iSly  +  S2)  y  } 


It  follows  that  the  ijth  term  of  {a^jy}{a^^^}^can  be  given  by 


termiij)  =  AY (ij)  ^(Sly  +  S2)  ' 'y 

r  =  1 


(A.80) 


(A.81) 


where  AY(iJ)^  mY(iJ)^  and  nY(iJX  are  given  in  Table  A.l  (LIV90).  Upon  substitution  of 
Eq.  A.81  into  Eq.  A.79,  the  stiffness  matrix  term  in  the  ith  row  andyth  column  for  theyi'th 
spar  may  be  expressed  as 

17/  Nh  Nh  5 


js  ~  4  V  L  ) 

ik~l  jk-l  r=l 


^spOj^ 


Jjy==^yL 


+A 


VA  syjf  -  sy^JJy^^yt 

L  \  fy=% 


{Sly  +  S2)  y  dy 


(A.82) 


syj^  -  syj^)h=^yL 


)C 


(51y  +  52)  y 


Table  A.l:  Coefficients  and  Powers  for  Curvature  Along  a  Spar  Line 


D 

AYUJX 

mYiiJX 

nY{ijX 

m 

imqi)imqj)(mqrl)(mqj-l)Sl^ 

mqrA  +  mqj 

nq,  +  nq, 

2 

2{mqi){mqj){mqr  1  ){nqj)S  1  ^ 

+  2{mqi){mqj)imqfl)(nqi)S\^ 

mqi  +  mqj-3 

nq^nqj-l 

3 

{mq^(nq^  {mqr  1  )(nqr  1  )6'  1  ^ 

+  4{mqi){mqj){nqd{nqj)S\^ 

+  {nq;)(mqj){nqi-  l){mqf  1)51 

mqf2  +  mqj 

nqi  +  nqf2 

4 

2imqi)(nqi)inqj- 1  )(nqj)S  1 
+  2inqi)(mqj)(nqrl)inqj)Sl 

mqi+mqj-l 

nq,  +  nqj-3 

5 

(nqi)(nqj)inqr  i){nqj- 1 ) 

nuji  +  mqj 

nq,  +  nqj-4 
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where 


m  =  mh-/^  +  mhji^ 

n  =  nY(iJ)  ^  +  nh.j^+nhji^ 


(A.83) 

(A.84) 


Noting  that  Eq.  A.82  is  composed  of  linear  combinations  of  the  integral  family  Isp{m,n) 
defined  in  Eq.  A.70,  the  final  expression  for  the  ijth  term  of  is 

Nh  Nh  5 

ZJ  z, 

"jk. 


r-l 


sp 


ik=^ljk=l  r-l  (A,85) 

^^spOj/sP  ^splj/sp  ^  +  1)  ] 


(LIV90).  The  stiffness  matrix  for  all  the  spars  combined  [K^p\toi  is  obtained  by  summing 
the  contributions  of  all  the  spars  given  by  Eq.  A.  85.  Because  the  wing  is  assumed  to  be 
symmetiic  with  respect  to  the  x-y  plane,  the  total  spar  contribution  can  be  found  by  sum¬ 
ming  the  contributions  from  the  spars  on  the  upper  wing  surface  and  then  multiplying  by 
2.  This  matrix  is  then  appropriately  merged  into  the  global  stiffness  matrix 


A.9  Mass  Matrix  Contributions  of  a  Rib 


Neglecting  rotary  inertia  the  contribution  to  the  kinetic  energy  T  of  an  element  of 
length  dx  of  the  yrth  rib  is 


(A.86) 


where  p^y  is  the  constant  material  density  of  the  rib,  and  the  cap  area  is  a  linear  func¬ 
tion  of  jc  (L1V90).  The  rib  cap  area  defined  in  Eq.  1.4  is  given  by 


+  A 


(A.87) 


Using  Eq.  A.8  for  wix,y,t),  Eq.  A.86  becomes 
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{rxpR,ry^ 

(rxAorjft) 


Figure  A.4:  Rib  Geometry 


dTj^  =  -Pj/i  (x)  {qj  {aj{ aj  { qj  dx 


Integrating  Eq.  A. 88  over  the  length  of  the  rib  gives 

1 
2 


where  the  limits  of  integration  rxp  and  vxa,  as  shown  in  Figure  A.4,  are  given  by 


rXp  -  Fly +  F2 


rx.  -  A  ly^j/D  +  A2 


where 


FI 


rxpg-rxp^ 

^yR-^yi 


F2  = 


rXpi/yR-rXppryL 

n^R-ryi 


Al 


f^AR-^^AL 


(AM) 

(A.89) 

(A.90) 

(A.91) 

(A.92) 

(A.93) 


^yR-R^L 


(A.94) 
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A2  = 


^^AL^R-^^AR^L 


ryR-^yi 


(A.95) 


Equivalent  to  Eq.  A.89  is  the  expression 
1  2 


Tjr  =  2‘"  Pir!ZZ‘A<^’‘'>rkA‘lJ  {«»>{«»}  fO* 


(A.96) 


where  to  is  the  natural  vibrational  frequency  of  the  rib.  Application  here  of  the  variational 
extremum  condition  from  Eq.  1.10  yields  the  mass  matrix  for  theyrth  rib  given  as 

{a,}  {aydx 


(A.97) 


Substituting  Eq.  A.87  for  A  ,  gives 


(A.98) 


Since  y  is  fixed  at  a  constant  value  y^yg  along  the  rib,  the  mass  matrix  term  in  the  ith  row 
and  yth  column  for  the  jrth  rib  may  be  expressed  as 


[Mrb  (*>y)  ]^v  -  Pjr[^rhOjyRlBlx=rxy’^^^ 


(A.99) 


where 


m  =  mq.  +  mq. 
n  -  nq.  +  nqj 

Let  us  define  a  family  of  integrals  where 

Then  the  final  expression  for  the  iyth  term  of  [M^^ljr  is 

i^rb  1  ,v  ~  Pjr  ^^rbQ.  ^ RB  ’^^rbljRB  +  1’  1 

J'  ''  Jr  jr 


(A.100) 

(A.101) 


(A.102) 


(A.103) 


(LIV90).  Appendix  C  discusses  how  the  family  of  integrals  is  analytically 
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derived.  The  mass  matrix  for  all  the  ribs  combined  is  obtained  by  summing  the 

contributions  of  all  the  ribs  given  by  Eq.  A.  103.  Because  the  wing  is  assumed  to  be  sym¬ 
metric  with  respect  to  the  x-y  plane,  the  total  rib  contribution  can  be  found  by  summing 
the  contributions  from  the  ribs  on  the  upper  wing  surface  and  then  multiplying  by  2.  This 
matrix  is  then  appropriately  merged  into  the  global  mass  matrix  [Mgi^i,]. 


A.10  Stiffness  Matrix  Contributions  of  a  Rib 

The  contribution  to  the  potential  or  strain  energy  [/  of  an  element  of  length  dx  of  the 
jrth  rib  is 


dUj^  = 


(A.  104) 


where  Ej^  is  the  longitudinal  modulus  of  elasticity,  the  cap  area  is  a  linear  function  of 
X,  and  h(x)  is  the  wing  depth  distribution  along  the  rib  line  (LIV90).  The  depth  distribu¬ 
tion  may  be  expressed  as 


h{x)  -  ^  H.j^  ■  X 


(A.  105) 


since  y  is  equal  to  a  constant  y^jg  along  the  length  of  the  rib.  Using  Eq.  A.8  for  w(x,y,t), 
Eq.  A.  104  becomes 


(A.  106) 


Integrating  this  over  the  length  of  the  rib  gives 


2  T 


(A.  107) 


where  the  limits  of  integration  have  been  defined  in  Eqs.  A.90  and  A.91.  Application  here 
of  the  variational  extremum  condition  from  Eq.  1.10  yields  the  stiffness  matrix  for  theyrth 
rib  given  as 


(A.108) 
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Substituting  Eq.  A.87  for  A  .  and  Eq.  A.  105  for  h{x)  here  gives 


=  Sr 


F 


Nh  Nh 

Z 

ik=l  jk=l 

li  II 

•  (A 

yRiB 


(A.  109) 


•  (^rhOj+^rhxX)  dx 


V£)l 


Now  differentiating  Eq.  A.  9  twice  with  respect  to  x  and  setting  y  equal  to  a  constant  ygig 
gives 


{a 


W^XX 


T 

}  = 


mqi-2  nq^ 


(A.  110) 


Upon  substitution  of  Eq.  A.  110  into  Eq.  A.  109  the  stiffness  matrix  term  in  the  /throw  and 
7th  column  for  the  jrth  rib  may  be  expressed  as 

^  Nh  Nh 

ik=ijk=i  (A.  Ill) 


rbO 


y  rib]  x=rXr, 


•x=rx.  m  ,  . 

dx  +  A 


rhl 


n  rx=rx^ 

•  ^  RIB  ]x=rXp 


m  +  1 


where 


AX(iJ)  =  mq.{mq-\)mq.{mq.-\) 

^  ^  (A.112) 

m  =  mq-  +  mq.-A  +  mh;.+  mh-, 

'  ^  (A.  113) 

n  =  nq.  +  nq-  +  n/i-  +  nh.. 

'  ^  (A.  114) 

Noting  that  Eq.  A.lll  is  composed  of  linear  combinations  of  the  integral  family  l^B(m,n) 
defined  in  Eq.  A.  102,  the  final  expression  for  the  ijth  term  of  [K^bljr  is 

^  Nh  Nh 

ik=l  jk=\ 

^^rbOj/RB  '^^rblj/RB  +  1’  ^)  ] 


(A.  11 5) 
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(LIV90).  The  stiffness  matrix  for  all  the  ribs  combined  is  obtained  by  summing  the 
contributions  of  all  the  ribs  given  by  Eq.  A.lll.  Because  the  wing  is  assumed  to  be  sym¬ 
metric  with  respect  to  the  x-y  plane,  the  total  rib  contribution  can  be  found  by  summing 
the  contribution  from  the  ribs  on  the  upper  surface  and  then  multiplying  by  2.  This  matrix 
is  then  appropriately  merged  into  the  global  stiffness  matrix  [Kgig^. 


Appendix  B: 

Analytical  Foundations  of  an  Equivalent 
Plate  Model  Using  FSDPT 


B.1  Overview 

In  this  appendix  are  derived  the  equations  necessary  to  construct  mass  and  stiffness 
matrices  for  wing  structural  components  based  on  First  Order  Shear  Deformation  Plate 
Theory.  This  appendix  starts  with  a  discussion  of  the  assumptions  of  the  theory  and  the 
ensuing  strain-displacement  relationships.  Displacement  functions  are  then  defined.  Fol¬ 
lowing  this  foundation  are  the  derivations  of  the  mass  and  stiffness  matrix  equations  for  a 
single  skin  layer.  These  are  followed  by  the  derivations  of  the  mass  and  stiffness  matrix 
equations  for  a  single  spar  and  a  single  rib.  Last  are  the  derivations  of  the  mass  and  stiff¬ 
ness  matrix  equations  for  a  single  spar  web  and  a  single  rib  web. 


B.2  Assumptions 

First  Order  Shear  Deformation  Plate  Theory  (FSDPT)  differs  from  Classical  Plate 
Theory  in  that  the  Kirchoff-Love  assumption  is  not  employed.  Rather,  it  is  assumed  that 
plane  sections  normal  to  the  plate’s  midsurface  remain  plane  but  not  necessarily  normal  to 
that  surface  after  deformation.  This  is  analogous  to  Timoshenko  beam  theory  (RE84). 
Hence,  the  shear  strains  and  may  not  be  ignored.  It  is  assumed  that  the  out  of  plane 
displacements  are  small,  and  there  may  be  in-plane  displacement  of  points  at  the  plate’s 
mid-surface.  It  is  assumed  that  the  transverse  normal  strain  is  negligible  since  the 
transverse  deflection  does  not  vary  through  the  thickness  of  the  plate.  This  requires  that 
the  transverse  normal  displacement  w  be  a  function  of  x  and  y  only. 

There  is  no  requirement  that  the  wing  structure  must  be  symmetric.  Spar  webs  and 
rib  webs  are  included  in  the  model  along  with  the  skins,  spars  caps,  and  rib  caps.  It  is  not 
assumed  in  the  general  case  that  skin  thickness  is  small  compared  to  the  wing  depth. 
When  this  is  assumed  then  it  may  also  be  assumed  that  all  the  skin  layers  in  each  of  the 
upper  and  lower  skins  have  the  same  depth  distribution  (LIV194). 
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B.3  Strain-Displacement  Relationships 

The  established  assumptions  allow  the  displacements  in  the  x,  y,  and  z  directions 
respectively  to  be  approximated  by 


u  (x,  y,  z)  =  Uq  (x,  y)  +  z\\f^  (x,  y) 
V  (x,  y,  z)  =  Vq  (x,  y)  +  zXKj,  (x,  y) 
w{x,y,z)  =  >VQ(A:,y) 


(B.l) 

(B.2) 

(B.3) 


where  Uq,  Vq,  Wq  are  the  x,y,z  displacement  components  of  a  point  along  the  reference  mid¬ 
surface,  and  \14  and  \|/^  are  rotations  of  a  line  element,  originally  perpendicular  to  the  mid¬ 
surface  plane,  about  the  y  and  x  axes  respectively  (RE84,LrV194).  The  associated  strains 
are  given  by 


du 


dx  dx  ^dx 


(B.4) 


^3')'  dy  dy 


(B.5) 


_  du  ^v  _  ^^0 

dy  dx  dy  dx 


dy  dx 


(B.6) 


dVUr 


(B.7) 


dWr 

y  =  v  +— 

‘y^  dy 


(B.8) 


(LIV194,RE84).  These  displacement  and  strain  relationships  form  the  basis  for  the  mass 
and  stiffness  matiix  formulations. 
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B.4  Displacement  Functions 

To  use  the  Ritz  method  described  in  Chapter  1,  it  is  necessary  to  approximate  each  of 
the  5  x,y,t  dependent  deformation  fields  by  a  polynomial  series  given  in  vector  form  by 

T 

uAx,y,t)  =  {aj(x,y)} 


vAx,y,t)  =  {a^(x,y)} 


^^ix,y,t)  =  {a^{x,y)}  {q^it)} 

T 

\\f  (x,y,t)  =  {a4(x,y)}  {^4(0} 


WQ(x,y,t)  =  {ag(x,y)}  {q^it)} 


(B.IO) 


(B.ll) 


(B.12) 


(B.13) 


where 


T  mql;  nq\: 

{a^{x,y)}  =  {,..,x  y  , ...} 


(B.14) 


and  similar  expressions  are  used  for  {a2V,  [^4}^,  and  {a5}^(LrV194).  The  column 

vectors  {i?2}.  {^3}.  {^4}>  and  {^5}  contain  the  polynomial  coefficients  which  are  the 
generalized  displacements  solved  for  in  the  Ritz  formulation.  These  polynomials  have 
Nqi,  Nq2,  Nq-^,  Nq4,  and  Nq^  terms  respectively,  and  the  x  and  y  powers  are  predetermined. 


B.5  Mass  Matrix  Contributions  of  a  Skin  Layer 

Let  the  vectors  {^1},  {^2}^  [q^)^  and  {^5}  be  combined  into  one  column  vector 

of  generalized  displacements  {  ^}  where 


{q]  =  {9p^2’^3’?4’95} 


(B.15) 


Substituting  Eqs.  B.9-B.13  into  Eqs.  B.1-B.3  allows  the  m,v,w  components  of  deformation 
to  be  written  in  terms  of  {^}.  This  is  given  by 
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u 

V 

w 


{ +^[-5^1  1  {(lit) } 


(B.16) 


.1  \  r  mSO{r,i)  nSO{r/i) 

where  Sq  and  are  matrices  containing  polynomial  terms  or  the  lorm  x  y 

and  ('■’ ('■’')  respectively.  Partitioned  into  subvectors  of  dimension  1  x  Nq, 
(5'=1,2,..,5),  the  matrices  are  given  as 


0  0  0  0 
0  4  0  0  0 
0  0  0  0  al 


(B.17) 


[S,  ix,y)] 


0  0  ^3  0  0 

0  0  0  aj  0 

0  0  0  0  0_ 


(B.18) 


where  both  are  3  x  in  dimension  with 


^‘Itot  ~  Nq^+  Nq2  +  Nq^  +  Nq^  +  Nq^ 


(B.19) 


based  on  the  length  of  the  subvectors. 

Now  the  contribution  to  the  kinetic  energy  T  of  an  infinitesimal  skin  element  dx  by  dy 
of  the  jlth  skin  layer  is 


ii 

T 

u 

V 

.  w  > 

>  < 

V 

.  w  . 

f  dxdy 


(B.20) 


where  p^,  is  the  constant  material  density  of  the  skin  layer,  and  tji{x,y)  is  the  thickness  of  the 
skin  layer  given  in  Eq.  1.2  as 


tji  (x,  y) 


‘'"■ji 

I?}-. 


X  y 


(B.21) 
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(LIV194).  Substituting  Eq.  B.16  into  Eq.  B.20  gives 


dTji  =  ^^Pjitji(x,y){qf[slSQ  +  z[sls^)  +  z[s%j  +  z\s\s^j]{q}dxdy 


(B.22) 


Integrating  this  over  the  whole  trapezoidal  panel  area  yields 


(B.23) 


Equivalently, 


Tji  =  \^'^9ji\\tjiix,y)  + +  {q}dxdy 


(B.24) 


where  co  is  the  natural  vibrational  frequency  of  the  skin  layer.  The  variational  extremum 
condition  from  Eq.  1.10  applied  here  becomes 


Oi  ;/ 

=  0 
dq; 


(B.25) 


yielding  the  mass  matrix  for  the  y/th  skin  layer  given  as 


J  jl  =  Py7 j  j ^jl  \  ^0^1  j  +  J  +  ^  V *^[^1  j  ] 


(B.26) 


(LIV194).  Substitution  of  Eq.  B.21  for  tji{x,y)  gives 


^^skhi  ~  '‘y 


k=i  yx 


(B.27) 


The  variable  z  is  assigned  the  value  of  the  depth  distribution  of  they/th  skin  layer  given  in 
Eq.  1.1  by 

Nhj, 

^  mh-t  nh,u 

hjiix,y)  =  y 


(B.28) 
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ml  =  mt,+mSO{r,i)  +mSO(r,j) 

(B.31) 

nl  =  nt,+nSO{r,i)  +nSO(r,j) 

^  (B.32) 

m2  =  mt.+mh.,+mSO{r,i)  +  mSl  (r,j) 

^  ''  (B.33) 

n2  =  nt,+nh;.+nSO(r,i)  +nSl(r,j) 

*  (B.34) 

m3  =  mt.  +  mh..  +mSl  (r,  i)  +  mSO  (r,j) 

*  (B.35) 

n3  =  nt,+nh,.+nSl{r,i)  +nSO(r,j) 

*  (B.36) 
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m4  =  mt/^  +  +  mS  1  ( r,  i)  +  mS  1  (r,  j) 


(B.37) 


n4  -  ntj^  +  nh-i^  +  nhji^  +  nSl  (r,  i)  +nSl  (r,j) 


(B.38) 


Noting  that  Eq.  B.30  is  composed  of  linear  combinations  of  the  integral  family  IjR{m,n) 
defined  in  Eq.  A.21  and  discussed  in  Appendix  C,  the  final  expression  for  the  iV’th  term  of 
is 

Nt.,  3 

k=\  r=l 

Nh., 

+  \ H.,  (Ij.„(m2,n2)  + 1 {m3,n3)) 

(B.39) 

Nh.,Nhj, 
ik=l  jk=l 

The  mass  matrix  for  all  the  skin  layers  combined  is  obtained  by  summing  the  con¬ 

tributions  of  all  the  layers  given  by  Eq.  B.39.  This  matrix  is  then  appropriately  merged 
into  the  global  mass  matrix 


B.6  Stiffness  Matrix  Contributions  of  a  Skin  Layer 


Each  skin  layer  is  treated  as  a  plane  stress  panel  for  which  the  only  strains  of  concern 
are  Byy,  and  Substituting  Eqs.  B.9-B.13  into  Eqs.  B.4-B.6  allows  these  strains  to  be 
written  in  terms  of  the  generalized  displacements  {q)  defined  in  Eq.  B.15.  This  is  given 
by 


{ [/?oU,y)]  +z[^iU.y)]}  {qit)} 


(B.40) 


where  Rq  and  are  matrices  containing  polynomial  terms  of  the  form 
C{r,i)x  y  and  C{r,i)x  y  respectively.  Partitioned  into 

subvectors  of  dimension  1  x  Nq^  (5’=1,2,...,5),  these  matrices  are  given  by 
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[Ro(x,y)]  = 


T 

^I’x 

0 

0  0 

0 

T 

0  0 

T 

T 

0  0 

0  0 

T 

^3  ’jc 

0 

0  0 

0 

T 

^4’y 

0  0 

T 

^3  ’y 

T 

^4’x 

[R^ix,y)]  = 


where  the  notation  for  the  partial  derivatives  is  defined  by  the  examples 
T 


a 


1  ’X 


a  ,  ,  ^  r  ,  mqli-l  nqli 

=  {...,mq\^-x  y  , ...} 


(B.41) 


(B.42) 


(B.43) 


and 


T  a  r  1  ^  r  1  n?l,-l 

=  {-^nq\.^-x  y  ,...} 


(B.44) 


Both  matrices  are  3  x  Nq^^,  in  dimension  with  Nq^^,  having  been  defined  in  Eq.  B.  19. 

Now  the  contribution  to  the  potential  or  strain  energy  t/  of  a  infinitesimal  skin  ele¬ 
ment  dx  by  dy  of  the  jlth  skin  layer  is 


T 

r  ^ 

e 

e 

XX 

XX 

^  ^yy 

^yy 

Y.V  J 

f  dxdy 


(B.45) 


where  [Q]ji  is  the  material  constitutive  matrix  of  the  skin  layer,  3  x  3  in  dimension,  as 
defined  in  Eqs.  A.26-A.37,  and  tji(x,y)  is  the  thickness  of  the  skin  layer  given  by  Eq.  B.21 
(LIV194).  Substituting  Eq.  B.40  into  Eq.  B.45  gives 


{g'idjcdy 


(B.46) 
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Integrating  this  over  the  whole  trapezoidal  panel  area  yields 


The  variational  extremum  condition  of  Eq.  1.10  applied  here  becomes 

=  0 

dq. 

yielding  the  stiffness  matrix  for  the  ;7th  skin  layer  given  as 

yx 

+z(  «[  [  ei  ;,«o  )  +  21 1 )  ]  *''>’ 

Substitution  of  Eq.  B.21  for  tji(x,y)  and  Eq.  B.28  for  the  variable  z  gives 


*»  . 

=  il 

k=l  ^yx 


mti  ntj  T  1 

X  V  [RolQ^jiRofxdy 


Nhj, 

ik=l  yx 


mt^  +  mh^,^  nt^  +  nh^ 


[[Rl[Q]jiR^  +  [R\[Q]jiRoj]^^dy 


Nhj,Nhj, 


+ 


u  Tj  C(^h  +  mhi;,  +  mhj,nt^+nh,,  +  nhJT  ,,  , 

[RtlQ]j,Kqdxdy 


ik=i  jk=l 


yx 


Let  the  r,/th  term  of  the  matiices  Rq  and  Ri  be 


(B.47) 


(B.48) 


(B.49) 


(B.50) 


.  h  r  ■\  nR0(r,i) 

RQ(r,i)  =  Roir,i)x  y 


(B.51) 
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(B.52) 

respectively.  Upon  substitution  of  Eqs.  B.51  and  B.52  into  Eq.  B.50,  the  stiffness  matrix 
term  in  the  ith  row  and  the  yth  column  for  the  jlth  skin  layer  may  be  expressed  as 
Ntji  3  3 

=  '^'^^fji[RQ(r,i)RQ(s,j)jjx^^y''^dxdy 

fc=ln=ly=l  L  yx 

Nhj, 

+  X  0^1  is,j)  jjx”'^y'‘^dxdy 

iic=i  ‘  (B.53) 

Nh., 

X  \\x”'^y''^dxdy 

ik=\  yx 

Nhj,Nhj, 

+  X  X  Wx^^^'^dxdy  1 

ik^\jk=\  ^  yx  -I 

where 


ml  =  mti^  +  mRQ  (r,  i)  +  mRO  (s,j) 
nl  =  nti^  +  nRO  (r,  i)  +  nRO  (s,j) 
m2  =  mtf^  +  mh^i^  +  mRO  (r,  i)  +  mR  1  is,j) 
n2  =  ntj^  +  nh.i^  +  nRO  {r,  i)  +nRl(s,j) 
m3  -  mtj^  +  mh-i^  +  mR  1  (r,  i)  +  mRO  (sj) 
n3  =  nti^  + nh^j^  +  nRl  (r,i)  +nR0{s,j) 
m4  =  mtj^  +  rnh^!^  +  mhjj^  +  mR  1  (r,  i)  +  mR  1  {s,j) 
nA  =  ntj^  +  nh-j^  +  nhjj^  +  nR  \  {r,  i)  +nR  \  (sj) 


(B.54) 

(B.55) 

(B.56) 

(B.57) 

(B.58) 

(B.59) 

(B.60) 

(B.61) 
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Noting  that  Eq.  B.53  is  composed  of  linear  combinations  of  the  integral  family  lTn(m,n) 
defined  in  Eq.  A.21,  the  final  expression  for  the  ijth  term  of  is 

Ntj,  3  3 

Tji  [~^Q  (r,  i)  Rq  is,j)  Ij-if  (fttl,  ra  1) 

k=l r-1 s=l 

Nhj, 

ik=l  (B.62) 

Nhji 

+  (^’  ^’)  ^0  (■^’2)  ("*3,  n3) 

/A:=l 

Nh,,  Nh., 

+  S  S  ^7-/?  ("^4,  n4)  ] 

/A:=l  yA:=l 

The  stiffness  matrix  for  all  the  skin  layers  combined  is  obtained  by  summing  the 

contributions  of  all  the  layers  given  by  Eq.  B.62.  This  matrix  is  then  appropriately  merged 
into  the  global  stiffness  matrix 

B.7  Mass  Matrix  Contributions  of  a  Spar 

The  contribution  to  the  kinetic  energy  T  of  an  element  of  length  dx\  of  the  jsth  spar  is 

T 


where  Tj  is  the  coordinate  along  the  spar  axis  rotated  from  the  y  axis  by  an  angle  A,  is 
the  constant  material  density  of  the  spar,  the  cap  area  A  is  a  linear  function  of  T|,  and  the 

*js 

velocity  vector  is  the  time  derivative  of  the  displacement  vector  defined  in  Eq.  B.16. 
Referring  to  Figure  A.3,  aU  T)  dependence  of  Eq.  B.63  may  be  changed  to  y  dependence. 
Substituting  Eq.  A.55  for  dr\  into  Eq.  B.63  and  allowing  the  cap  area  to  be  expressed  as  a 
linear  function  of  y  gives 

'  T  ' 

u  u 

V  V 

w  )  I  |{; 


(B.64) 
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Substituting  Eq.  B.16  into  this  gives 


+z 


(B.65) 


Here  the  spar  line  equation  from  Eq.  A.57  may  be  employed  to  express  x  within  the  matri¬ 
ces  ^0  “ind  5'i  as  a  function  of  y.  Therefore  the  polynomial  terras  contained  in  Sq  and  take 
the  form  (5^1y-i-52)  y  and  (51y-i-52)  y  respectively. 

Integrating  Eq.  B.65  over  the  length  of  the  spar  now  gives 


T.  =  Ip.f - ^ - V  '^''A(y)^  + 

2^J\syj^-syJjy=^yL  L  o  o  v  o  i 


Equivalently, 


(B.66) 


(B.67) 


where  co  is  the  natural  vibrational  frequency  of  the  spar.  Application  here  of  the  varia¬ 
tional  extremum  condition  from  Eq.  1.10  yields  the  mass  matrix  for  the  y.yth  spar  given  by 


js  syjf  -  syi^Jiy=^yL 


+4'SX)  +  z'(5[s,)]rfy 


(B.68) 


Substituting  Eq.  A.56  for  the  cap  area  and  Eq.  A.74,  the  spar  depth  distribution  in  terms  of 
y,  for  the  variable  z  gives 
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js 


ik=l 


(B.69) 


ik~\  jk=\ 


The  mass  matrix  temi  in  the  ith  row  and yth  column  for  theysth  spar  may  be  expressed  as 

3 


=  P„ 


ml  n\ 


:]s  K:(^.po, +4,,/)  -«)■■■  v-i. 


*  E  »JZ[  ^  ” V"  ^  (Sly  *S1^ " dy 

ik=\ 

- 1 E  A.n:  (^..0,.  -  ] 


ik-\  jk=l 


(B.70) 


where 


ml  =  mSO(r,i)  +mSO(r,j) 
nl  =  nSO{r,i)  +nSO{r,j) 
m2  =  mh^i^  +  mSO{r,i)  +mS\  (rj) 
n2  =  nh-i^  +  nSO  (r,i)  +nSlir,j) 
m3  =  rnh-j^  +  mSl  (r,  i)  +  mSO  {r,j) 
n3  =  nh-j^+nSl  (r,i)  +nSQ(r,j) 
m4  =  mh-i^  +  mhjj^  +  mS\  (r,  i)  +  mS  1  (r,  j) 


(B.71) 

(B.72) 

(B.73) 

(B.74) 

(B.75) 

(B.76) 


(B.77) 
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n4  =  nh.i^  +  nh-j^  +  nSl  {r,i)  +nSl  (r,j) 


(B.78) 


Noting  that  Eq.  B.70  is  composed  of  linear  combinations  of  the  integral  family  Isp{m,n) 
defined  in  Eq.  A.70  and  discussed  in  Appendix  C,  the  final  expression  for  the  ijth  term  of 
[Msplji  is 

3 

r=l 

Nhj, 

ik~l 

(B.79) 

Nhj, 

+  X  ^js^/^spij,  (4/>  ("*2,  n2+l)+  I^p  (m3,  n3  +  1) ) 

ik=^l 

m„Nh„ 

+  I  D)  ] 

/A:=l  7^=1 


The  mass  matrix  for  all  the  spars  combined  is  obtained  by  summing  the  contribu¬ 

tions  of  all  the  spars  given  by  Eq.  B.79.  This  matrix  is  then  appropriately  merged  into  the 
global  mass  matrix 


B.8  Stiffness  Matrix  Contributions  of  a  Spar 


The  contribution  to  the  potential  or  strain  energy  f/  of  an  element  of  length  dx\  of  the 
75th  spar  is 


(B.80) 


where  is  the  longitudinal  modulus  of  elasticity,  the  cap  area  A  is  a  linear  function  of 

^rjs 

the  spar  line  coordinate  T],  and  is  the  normal  strain  along  the  spar  axis.  Referring  to 
Figure  A.3  and  using  standard  tensor  transformation  rules,  the  normal  strain  in  the  r\  direc¬ 
tion  may  be  written  in  terras  of  strains  in  the  x  and  y  directions  by 
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where 


Sti 


{TRf] 


XX 


yy 


*xy 


(B.81) 


{77?}^=  { sin^A,  cos^A,  sinAcosA}  =  {^^A,  c^A,  ^AcA) 


(B.82) 


The  spar  cap  area  may  easily  be  expressed  as  a  linear  function  of  y  as  shown  in  Eq.  A.56. 
Substituting  Eq.  A.55  for  dx]  and  Eq.  B.81  for  into  Eq.  B.80  gives 


— - — y^yhp 

2  ^\syj^-syJ 


yy 


T, 


{TR} {TR} 


xy 


"yy 


•  xy 


dy 


Define 


[Qsp]  =  {TR}  {TA}  = 


(s^A)  (^2Ac2A)  (^3AcA) 

(s^Ac^A)  (c'^A)  (^Ac^A) 

(s^AcA)  (^Ac^A)  {s^Ac^A\ 


(B.83) 


(B.84) 


Substituting  Eqs.  B.40  and  B.84  into  Eq.  B.83  gives 


[e.,1  ''o)  +  [e,pl R,) ]  {«> dy 


(B.85) 


Here  the  spar  line  equation  from  Eq.  A.57  may  be  employed  to  express  x  within  the  matri¬ 
ces  Rfj  and  Ai  as  a  function  of  y.  Then  the  polynomial  terms  contained  in  Rq  and  Ri  take  the 
form  C(r,i)  (Sly  +  S2)  y  and  C  (r,  i)  (Sly  +  S2)  y 

respectively.  Integrating  Eq.  B.85  over  the  length  of  the  spar  yields 
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U  -  -E  f  — 

~  2  J\syj^-syJh='^yL 


(,)  ^[rI  [e,,]  + <^0  [2,,]  /?! 


(B.86) 


+4  r\iQJ  Ro  J  +  z’l  R\iQsp^  RA']{q}dy 


Application  here  of  the  variational  extremum  condition  from  Eq.  1.10  yields  the  stiffness 
matrix  for  the  ysth  spar  given  as 


^^Ksyj^-sy^^ 


AR\{QJRA^z\R\[Q,p\RAyy 


(B.87) 


Substitution  of  Eq.  A.56  for  the  cap  area  and  Eq.  A.74  for  the  variable  z  gives 


^  ^  ji  V  “  ^yi 


(B.88) 


+  X  l’'  ^  ^2)  "‘'V'’"  -f  A,,,. ,.)  r\ [Q,^]  R„dy 


+  y  YH-  H-  r"''(Slj.  +  i2) 

Zrf  Z-  JSjJy=sy^_  '' 


mhi^  +  mhji^  nh,^  +  nhj^ 


ik=\  jk=\ 


Eqs.  B.51  and  B.52  are  used  to  define  the  r,/th  term  of  the  matrices  Rq  and  R^  respectively 
with  X  replaced  by  Eq.  A.57.  Upon  substitution  of  Eqs.  B.51  and  B.52  into  Eq.  B.88,  the 
.stiffness  matrix  term  in  the  ith  row  and  the  yth  column  for  theysth  spar  may  be  expressed 


\ro  0  Ro  fyZll  ^^spoj, +^.pi/) 


(B.89) 


+ 1  HjJoC-.  0  R,  is.jl  (Sly  +  52)  "V^dy 


A  2  Hj,k,  (r.  i)  R„  IsJ)  f;;:''  (A,,o,_  +  J)  (5 b  +  52)  "^/^dy 


Nhj^Nhj^ 


S««  «/./>  <'■•  (51y  +52)”V%] 


,fc=l  yfe=l 


where 

ml  =  mRO(r,i)  +  mRO  (s,j) 

(B.90) 

nl  =  nRO(r,i)  +nRO(s,j) 

(B.91) 

m2  =  mh;.+  mRO  (r,i)  +mRl  is,j) 

(B.92) 

n2  =  nh.,+nRO{r,i)  +nRl{s,j) 

(B.93) 

m3  =  mh.,+mRl  (r,  i)  +mi?0(5,7) 

(B.94) 

n3  =  + /i/?l  (r, /)  +nRO(s,j) 

(B.95) 

m4  =  m/i,.  +  mh;^  +  mRl  (r,  i)  +mRl  (sj) 

^  ^  (B.96) 

nA  =  nh;,+nh..+nR\{r,i)  +nR\{s,j) 

^  (B.97) 


Noting  that  Eq.  B.89  is  composed  of  linear  combinations  of  the  integral  family  Isp{m,n) 
defined  in  Eq.  A.70,  the  final  expression  for  the  /  Jth  term  of  is 

3  3 

r=l  5'=1 

[Ro{r,i)Ro{s,j)  iA^pQjsp{ml,nl)  +A^^j,/5p(ml, nl  +  1)) 
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+  Z  (r,  i)  Ri  isj)  (A^pOj/sp  ("^2,  n2)  +  (m2,  n2  +  1) ) 

ik=i  “  ''  ''  (B.98) 

(A^pQjgpim3,n3)  +  (m3, n3  +  1)) 

1^=1 

m„Nhi, 

+  S  ('4,pO„V('"'',n4)  +A,^,^_V(m4,n4  +  l))  ] 

iA;=l  7^=1 

The  stiffness  matrix  for  all  the  spars  combined  is  obtained  by  summing  the  contri¬ 
butions  of  all  the  spars  given  by  Eq.  B.98.  This  matrix  is  then  appropriately  merged  into 
the  global  stiffness  matrix 


B.9  Mass  Matrix  Contributions  of  a  Rib 

The  contribution  to  the  kinetic  energy  7  of  an  element  of  length  dx  of  the  jrth  rib  is 


1 


T 

' 

u 

u 

<  “ 

>  < 

V 

V 

.  w  > 

.  w  . 

dx 


(B.99) 


where  pj^  is  the  constant  material  density  of  the  rib,  the  cap  area  A^^  is  a  linear  function 
of  X,  and  the  velocity  vector  is  the  time  derivative  of  the  displacement  vector  defined  in 
Eq.  B.16.  Substituting  Eq.  B.16  into  Eq.  B.99  gives 

dTjr  =  -jpj^Aix)  4- j  -1- z 

(B.lOO) 


Since  y  is  equal  to  a  constant  y^jB  over  the  length  of  the  rib,  the  polynomial  terms  con- 

1  r*  m50(r, /)  nSO{r,i)  ,  mS\{r,  i)  (r,  i)  .  , 

tained  in  and  take  the  form  x  and  x  y^jg  respectively. 

Integrating  Eq.  B.lOO  over  the  length  of  the  rib  now  gives 


(B.lOl) 
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where  the  limits  of  integration  rxp  and  rx^  are  defined  in  Eqs.  A.90  and  A.91  referencing 
Figure  A.4.  Equivalent  to  Eq.  B.lOl  is  the  expression 


(B.102) 


where  co  is  the  natural  vibrational  frequency  of  the  rib.  Application  here  of  the  variational 
extremum  condition  from  Eq.  1.10  yields  the  mass  matrix  for  theyrth  rib  given  as 


(B.103) 


Substituting  Eq.  A.87  for  the  rib  cap  area  and  Eq.  A.  105  for  the  variable  z  gives 


[M,,]  = p, 


Nhj^  /  A  /  A 

ik^i  (B.104) 

Nhj^Nhj^  /  N 


ik=l  jk=l 


The  mass  matrix  term  in  the  fth  row  and yth  column  for  the yrth  rib  may  be  expressed  as 


Wrb  ( h  j)  ]  -  PyvX  [^i?/sjx=r4  ^ 


n2  c^=rxA  .  .  .  .  ml  , 

(''rW,  ^  * 


n3  rx~rx.  .  .  .  .  m3  , 


Nhj^Nhj 


+2 


ik=l  jk=l 


(B.105) 
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where  the  x  and  y  powers  are  defined  in  Eqs.  B.71-B.78.  Noting  that  Eq.  B.105  is  com¬ 
posed  of  linear  combinations  of  the  integral  family  Igg(m,n)  defined  in  Eq.  A.  102  and  dis¬ 
cussed  in  Appendix  C,  the  final  expression  for  the  ijth  term  of  is 

3 

r=l 

Nhj, 

'^^^jr.^rho.  +/^g(m3,  n3)) 

<vfc=i  “  ''  (B.106) 

Nh^r 

+  ^  '''  ^RB  ^ 

ik=\ 

Nhj^Nhj^ 

+  X  +A^^,,^(m4+l,n4))  ] 

ik=ljk=\ 


The  mass  matrix  for  all  the  ribs  combined  is  obtained  by  summing  the  contribu¬ 

tions  of  all  the  ribs  given  by  Eq.  B.106.  This  matrix  is  then  appropriately  merged  into  the 
global  mass  matrix 


B.10  Stiffness  Matrix  Contributions  of  a  Rib 


The  contribution  to  the  potential  or  strain  energy  [/  of  an  element  of  length  dx  of  the 
jrth  rib  is 


du.^  =  2  V 


,  6^  dx 

rbj^  XX 


(B.107) 


where  is  the  longitudinal  modulus  of  elasticity,  and  the  cap  area  is  a  linear  func¬ 
tion  of  X  Using  Eqs.  B.4,  B.9,  and  B.  1 1  let  us  define 

+z[FiUy)]} 

(B.108) 


where 


T  T  T 

W  = 


[yoU^)]  =  [a[,,o] 
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(B.109) 

(B.llO) 

(B.lll) 


Both  matrices  To  and  Fj  are  partitioned  into  2  subvectors,  1  x  Nqi  andl  x  Nq^  in  dimension 
respectively.  Since  y  is  equal  to  a  constant  along  the  length  of  the  rib,  and  Fj  con- 

myo.  fiYO.  mKl.  nYl-  .  .  1/^,0 

tain  polynomial  terms  of  the  form  C^x  C^x  .  Substitutmg  Eq.  B.108 

into  B.  107  gives 

dUj,  =  {«)  +=:(  'Ti'o)  +  ffj'i)]  «> 

(B.112) 


Integrating  this  over  the  length  of  the  rib  gives 

(B.113) 

where  the  limits  of  integration  have  been  defined  in  Eqs.  A.90  and  A.91.  Application  here 
of  the  variational  extremum  condition  from  Eq.  1.10  yields  the  stiffness  matrix  for  theyrth 
rib  given  as 


(B.114) 


Substituting  Eq.  A.87  for  the  rib  cap  area  and  Eq.A.105,  the  rib  depth  distribution,  for  the 
variable  z  gives 


+  S  [(  >^>'l )  +  (  if  '’o)] 


ik=\ 

Nhj^Nhj^ 


(B.115) 


i^=l  jk=\ 


•x=rx.  mhii^  +  mhji^  nh^^  +  nh 
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Let  the  rth  term  of  the  single  row  matrices  Yq  and  Fj  be 

-  mm  tiYO, 

4  =  V  > 


^  mYL  nYli 


(B.116) 

(B.117) 


Upon  substitution  of  Eqs.  B.116  and  B.117  into  Eq.  B.115,  the  stiffness  matrix  term  in  the 
ith  row  and  jth  column  for  the  jrth  rib  may  be  expressed  as 


Nhjr 


+  X  i^bO^+^rbl/)  "'dx 

ik=l 
Nhjr 

+  X  "jr/iJo/lUT-Zl 


(B.118) 


ik=l 


Nhj^Nhj^ 


-  s  X  cn:  ] 


ik=l  jk=l 


where 


ml  =  mFO. +  myO^- 

(B.119) 

nl  =  nY0.  +  nYOj 

(B.120) 

m2  =  mh.j^  +  mY0-  +  mYlj 

(B.121) 

n2  =  nh^i^  +  nY0.  +  nYlj 

(B.122) 

m3  =  mh.i^  +  mYl  ■  +  mYOj 

(B.123) 

n3  =  nh-i^  +  nYl.  +  nYOj 

(B.124) 

124 


m4  =  mh..+mh-,^  +  mYl+mY\. 

Ik  jk  j  (B.125) 

«4  =  nh;.-\-nh..+nY\.  +  nY\. 

Ik  jk  j  (B.126) 

Noting  that  Eq.  B.  118  is  composed  of  linear  combinations  of  the  integral  family  IgB(m,n) 
defined  in  Eq.  A.  102  and  discussed  in  Appendix  C,  the  final  expression  for  the  ijth  term  of 
Myris 

[KrijiUj)]  Jr  -  ^rblj RB 

Nh.^ 

+  E  ^jrJoS  ^^rbojRB  «2)  +  (m2  +  1,  n2) ) 

ik=i  '  (B.127) 

Nhj, 

ik-1 

+  E  S«yr,Ar/iri,(^rM,/K«<'»4,n4)  +  (m4  +  1.  n4) )  ] 

/A:=l  jk=l 

The  stiffness  matrix  for  all  the  ribs  combined  is  obtained  by  summing  the  contribu¬ 

tions  of  all  the  ribs  given  by  Eq.  B.127.  This  matrix  is  then  appropriately  merged  into  the 
global  stiffness  matrix  [Kgi„^. 

B.11  Mass  Matrix  Contributions  of  a  Spar  Web  Layer 

A  spar  web  is  positioned  in  the  vertical  plane  between  parallel  spar  caps  on  the  upper 
and  lower  wing  surfaces.  This  plane  is  defined  to  be  the  T|-z  plane  where  T|  is  the  coordi¬ 
nate  along  the  corresponding  parallel  upper  and  lower  spar  axes  rotated  form  the  y  axis  by 
an  angle  A.  See  Figure  A.3  for  spar  line  geometry.  The  contribution  to  the  kinetic  energy 
T  of  an  infinitesimal  element  dx\  by  dz  of  the  y7th  spar  web  layer  is 

"  d^dz 

(B.128) 
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where  pji  is  the  constant  material  density  of  the  web  layer,  the  velocity  vector  is  the  time 
derivative  of  the  displacement  vector  defined  in  Eq.  B.16,  and  tji{T\,z)  is  the  thickness  of 
the  layer.  The  layer  thickness  may  be  expressed  as  a  linear  function  y  only.  Referring  to 
Figure  A.3,  all  T|  dependence  may  be  changed  to  y  dependence  by  substituting  Eq.  A.55 
for^iTl.  This  gives 


L 

^yR-^yi 


(B.129) 


where  the  parameters  L,  sy^,  and  syi  may  be  defined  from  either  the  upper  spar  cap  or  the 
lower  spar  cap,  and  the  layer  thickness  is  given  by 


(B.130) 


Substituting  Eq.  B.16  into  Eq.  B.129  gives 


dTi,  =  2?, 7 


^yR-^yi 


7/ 


+d  l  +  z 


(B.131) 


{q}dydz 


Here  the  spar  line  equation  from  Eq.  A.57  may  be  employed  to  express  x  within  the  matri¬ 
ces  So  and  5i  as  a  function  of  y.  Therefore  the  polynomial  terms  contained  in  Sq  and  take 
the  form  (51y  +  52)'”^°^''’'^y"^“^'‘’'^  and  (51y +  S2)'”^^^''’'^y"^^^'’’'^  respectively. 
Integrating  Eq.  B.130  over  the  area  of  the  web  layer  now  gives 


\_  (  L  n=sy^ 


syR-syi 


{q}dydz 


(B.132) 


where  the  limits  of  z  integration  ai'e  the  depth  distributions  of  the  lower  and  upper  spar 
caps  given  respectively  by 
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Nh 


u 


hfjUy)  =  y 


ik  -  1 

Nh, 


hi  ix,y)  =  ^  HTj^-x  y 


ik  =  1 


Equivalent  to  Eq.  B.  131  is  the  expression 

■  3-X,i^JC:C:S'-“ 

+z(5X)  +  z^(5[5j]  {^} 


(B.133) 


(B.134) 


(B.135) 


where  co  is  the  natural  vibrational  frequency  of  the  web  layer.  Application  here  of  the 
variational  extremum  condition  from  Eq.  1.10  yields  the  mass  matrix  for  they/th  spar  web 
layer  given  as 


jl 


9ji 


^yR-^yt 


sis. 


(B.136) 


The  z  integration  is  performed  first  analytically  using 


(hu(x,y) 


(B.137) 


Eq.  B.136  thus  becomes 


j I  ~  9ji 


^yR-^yi 


+\{hl  (x,  y)  -  h\  {x,  y) 


(B.138) 
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Substituting  Eq.  B.130  for  the  layer  thickness  and  Eqs.  B.133  and  B.134  for  the  spar  cap 
depth  distributions  gives 

rf^hj 


^^swh^ji  Py/ 


\-ik=\ 


t/p'=%  ,  ft  Yc^c  u, 


Nh.^  . 


ik=\ 

NK  Nh 


iyn^i\ny  x//  \ 


ik=l  jk-\ 


Nh^  Nh^ 


1 V  V  fy=^yx 
2  2-(  '■*  7*Jy=*yL 


mhL;j,  +  mhLjf^ 


f^=l  yA;=:l 

NhjjNhfjNh^ 


(B.139) 


U  U  U  P’=syg  mhU.^+mhUji^  +  mhU,^  nhU.i^  + nhUj^  + nhU,^^ 

3Zu  ^  ^^■^ik"jk^lk]y=sy^^  y 

ik=l  jk=l  lk=l 


1  v-'"  -,L  --L  fy=% 

-3  Z  Z  X  >■ 

■(^o„+v)(«[^,Y] 


7V/^^ 


/^=1  jk~l  lk=l 


Substituting  the  spar  line  equation  from  Eq.  A.57  for  all  x,  the  mass  matrix  term  in  the  ith 
row  and  yth  column  for  the  jlth  spar  web  layer  may  be  expressed  as 


Pji\ 


Nh, 


^yR-^yu% 


3  rNh^ 


^ik=l 


ik=\ 


NhyNhy 


(Sly.S2rV^^)dy 


ik=l  jk~l 


Nh^ 


ik=\  jk=l 

Nhy  NhyNhy 


(B.140) 


■is  s  s<.«xn:^,-vF‘> 


ik=l  jk=l  lk=l 
Nh^  Nh^  Nhj^ 

1  X"'  TJ^  TJ^  TJ^  0'=% 


1  jji^  jjL 

1  X  2.  2.  ^ik^jk^h 


A,  A,  ^  ^  ^  X  nL4  , 

7^0  +^iJ  (*^ly  +  '^2)  y  dy 


%v 


ik=l  jk=l  lk=l 


where 

mill  =  mhU .,  +  mSO  {r,  i)  +  mSO  (r,j) 

*  (B.141) 

nUl  =  nhU,.+nS0(r,i)  +nS0(r,j) 

'*  (B.142) 

mLl  =  mhL..+mS0(r,i)  +mS0{r,j) 

(B.143) 

nLl  =  nhL..+nSOir,i)  +nS0ir,j) 

""  (B.144) 

mU2  =  mhU-.+mhU..+mSO(r,i)  +  mSl  (r,j) 

Ik  jk  (B.145) 

nU2  =  nhU..+nhU,^  +  nSO{r,i)+nS\{r,j) 

^  (B.146) 

mL2  =  mhL;.+mhL..+mSO(r,i)  +mSl{r,j) 

tk  Jk  (B.147) 

nL2  =  nhL;.  +  nhL.,  +  nSO  {r,  i)  +nSl  (r,j) 

""  ^  (B.148) 

mU3  =  mhU;i,  +  mhU-^  +  mSl  (r,  i)  +  mSO  (rj) 

^  (B.149) 

nU3  =  nhUii^  +  nhUji^  +  nSl  ir,i)  +nS()ir,j) 


(B.150) 
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mL3  =  nthL-j^  +  mhLjj^  +  mS\  (r,  i)  +  mSO  (r,j) 


nL3  =  nhL-i^  +  nhLjj^  +  nSl  (r,  i)  +  nSO  (rj) 


mU4  =  mhU-i^  +  mhUj,^  +  mhUi/^  +  mSl  (r,  i)  +  mSl  (rj) 


nU4  =  nhU-i^  +  nhUji^  +  nhU^i^  +  nSl  {r,  i)  +nSl(r,j) 


mL4  =  mhL^i^  +  mhlj/^  +  mhL^i^  +  mSl{r,  i)  +mSl  (rJ) 


nL4  =  nhL.^j^  +  nhlj/^  +  nhLii^  +  nSl  (r,  i)  +  nSl  (rJ) 


(B.151) 

(B.152) 

(B.153) 

(B.154) 

(B.155) 

(B.156) 


Noting  that  Eq.  B.140  is  composed  of  linear  combinations  of  the  integral  family  Isp{m,n) 
defined  in  Eq.  A.70,  the  final  expression  for  the  ijth  term  of  is 


3  ^Nhi, 

Py/X 


r=l  \-ik=l 


"  , 

fojsp  (mUl,  n[/l  +  1)  J 


Nh^ 

■  X  "^1)  +  ^i/sp  + 1) 

ik-\ 


NhyNhy 


+3  X  X  Wsp  nU2)  +  imU2,  nU2+l) 


ik=l  jk=l 


Nhfj  Nhy 


ik=l  jk=l 


Nh,  Nh, 


4X  Vs? +  t,/sp(mL2,nL2+l)) 

ik=l  jk=l 


ik=i  jk-l 


(B.157) 
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Nh^jNhyNhu 

lk=l  jk=l  lk=l 
Nh^  Nh^  Nhj^ 

ik=\  jk=\  lk-\ 


The  mass  matrix  for  all  the  spar  webs  combined  is  obtained  by  summing  the  con¬ 

tributions  of  all  the  spar  webs  given  by  Eq.  B.157.  This  matrix  is  then  appropriately 
merged  into  the  global  mass  matrix 

B.12  Stiffness  Matrix  Contributions  of  a  Spar  Web  Layer 

A  spar  web  layer  in  the  T|-z  plane  is  treated  as  a  plane  stress  panel  where  the  only 
strains  of  importance  are  ,  and  .  From  the  assumptions  stated  in  section  B.2, 

^zz  neglected.  The  normal  strain  has  been  defined  in  terms  of  strains  in  the  x- 

y  plane  in  Eq.  B.81.  The  shear  strain  y  may  also  be  defined  in  terms  of  strains  in  the  x-y 
plane  using  standard  tensor  transformation  rules.  It  is  given  by 

V-nz  =  {sinA,  cosA}  1  =  {^'A,  cA}  \  [ 

I  Yj, J  I  ly,  >  (B.158) 


where  A  is  the  angle  of  rotation  of  the  q  axis  from  the  y  axis.  Combining  Eqs.  B.81,  B. 8 2, 
and  B.  158  we  may  write 


I  =  [TR^^]  {£} 

^TIZ 


(B.159) 


where 


{E}  -  {^xx'^yy’^xy’ 


(B.160) 


and  the  transformation  matrix  is  defined  by 
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[TR,A 


(s^A)  (c^A)  (5Ac'A)  0  0 

0  0  0  (sA)  (cA) 


(B.161) 


Using  Eqs.  B.4-B.13  and  Eq.  B.15,  Eq.  B.159  may  be  expressed  in  terms  of  the  general¬ 
ized  displacement  vector  {^}  in  the  form 


(B.162) 


where 


[W.U,y)]  = 


a\,^  0  0  0  0 

0  0  0  0 

«L  0  0  0 

0  0  ag  0 

0  0  0  al  al,y 


(B.163) 


[W,  ix,y)] 


0  0  0  0 

0  0  0  alyO 

0  0  4,^  4,^  0 

0  0  0  0  0 

0  0  0  0  0 


(B.164) 


The  matrices  Wq  and  Wj  contain  polynomial  terms  of  the  form  C(r,  i)  i)  ^nWQ{r,  i) 

and  C  (r,  i)  respectively  and  are  partitioned  into  subvectors  of  dimen¬ 

sion  1  X  Nq,  (5=1, 2,..., 5).  Both  matrices  are  5  x  Nq^^,  in  dimension  with  Nq,o,  having  been 
defined  in  Eq.  B.  19. 

Now  the  contribution  to  the  potential  or  strain  energy  (/  of  an  infinitesimal  element 
dr[  by  dz  of  the  jlth  spar  web  layer  is 


d^dz 


(B.165) 
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Figure  B.  1:  Positive  Rotation  of  Principal  Material  Axes  from  ^-z  Axes 


where  [Sly/  is  the  jlth  layer’s  constitutive  matrix,  and  tji(r\,z)  is  the  thickness  of  the  layer. 
The  constitutive  matrix,  as  discussed  in  Section  1.5,  has  the  form 


[Slj,  = 


Sii  Si6 

0l6  S66 


(B.166) 


where  its  components  are  determined  for  the  T|-z  axis  system  from  [Q]ji,  they/th  layer’s 
constitutive  matrix  referenced  to  its  own  principal  material  axes.  Assuming  they/th  layer 
to  be  orthotropic,  its  invariant  properties  may  be  used  to  determine  the  components  of 
[Sly/  when  the  layer’s  principal  axes  are  oriented  at  some  angle  P  to  the  Tj-z  axes  as 
shown  in  Figure  B.l  (J075).  The  components  of  [Sly/  may  be  written  as 


=  Cj  +  C2COS2P  +  C3COS4P 
S16  =  -^C2sin2p-C3sin4p 
See  =  C5-C3COS4P 


(B.l  67) 
(B.168) 
(B.169) 


for  which  the  invariants  Q,  C2,  Cj,  and  C5  are  defined  in  Eqs.  A.33-A.37  (J075).  The 
layer  thickness  may  be  expressed  as  a  lineai'  function  of  y  only  as  given  in  Eq.  B.l 30. 
Referring  to  Figure  A.3,  all  T|  dependence  of  Eq.  B.165  may  be  changed  to  y  dependence 
by  substituting  Eq.  A.55  for  dx].  This  gives 
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ir  L 


'n  -  2 


^yR-^yi 


tAy)\  [S]  ,7  \dydz 


(B.170) 


where  the  parameters  L,  syfi,  and  syi  may  be  defined  from  the  geometry  of  either  the  upper 
spar  cap  or  the  lower  spar  cap.  Substituting  Eq.  B.159  into  Eq.  B.170  gives 


If  L 


'n  -  2 


^yR-^yiy- 


tAy)  {^}dydz 


(B.171) 


Let  us  define 


[ej  =  (ra^l  ISljiLraj 


(B.172) 


where  [QJ  is  5  x  5  in  dimension.  Substituting  Eqs.  B.162  and  B.172  into  Eq.  B.171  gives 


(B.173) 


+4  [2,]  ^0  +  w[[Q^]wM{q}  dydz 


Here  the  spar  line  equation  from  Eq.  A.57  may  be  employed  to  express  x  within  the  matri¬ 
ces  Wq  and  Wi  as  a  function  of  y.  The  polynomial  terms  contained  in  Wq  and  Wj  take  the 
form  C(r,  j)  (‘S^ly  +  ‘S^2)  y  '  '  and  C(r,  0  (-S^ly  +  *S^2)  y 

respectively.  Integrating  Eq.  B.173  over  the  area  of  the  web  layer  now  gives 


(B.174) 


+4  [GJ  wJ  +  < [ GJ  wA]{q}  dzdy 


where  the  limits  of  z  integration  are  the  depth  distributions  of  the  lower  and  upper  spar 
caps  given  in  Eqs.  B.133  and  B.134.  Application  here  of  the  variational  extremum  condi¬ 
tion  from  Eq.  1.10  yields  the  stiffness  matrix  for  the  y/th  spar  web  layer  given  as 
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.  1  =(  ^  t 

L  swbiji  {syj^-  syJjy=syJf^L(x,y) 


(y)  K[<2,]Wo  +  z  <[!2J  Wj 


(B.175) 


+z  +  z[  W’[[(2j] 


Performing  the  z  integration  analytically  using  Eq.  B.137  gives 


i7 


(^. )')  -  h,  ix,  y) )  (  <  [e,] 

+!(/!?,  y)  -  y) )( ( <  [0,1  W', )  +  ( <  [aj  »o 

+\{hlix.y)  -hlu,y^][v/\{ QJ  W, J ] 7y 


(B.176) 


Substituting  Eq.  B.130  for  the  layer  thickness  and  Eqs.  B.133  and  B.134  for  the  spar  cap 
depth  distributions  gives 


['fw/.l 


yyR-syi, 


'Nh-y  .  . 

Lik^l 


=sy„  mhL-i^  nhLi^ 


y  “  T’o.  +  V 


0  ^^0 


]  W.  \dy 


NhyNhy 


U^Ury=syg  mhUi^  + mhUj^  nhU-i^  +  nhUj^ 


ik=\  jk=\ 


Nh,  Nh^ 


^0,  +  T,y  \W,[Q^]W,  U\  w;  ly,  dy 


1  mhLi^  +  mhLji^  nhLii^  +  nhLj^ 

■2I  >■ 

ik=l  jk~l 


K  +  ^1/  <  tej  W J  +  <  [0J  Wo  ]dy 


(B.177) 


NhyNhyNhy 


U  U  U  ry=syR  mhUi^^  +  mhUji^  +  mhUn  nhU.,^  +  nhUj^  +  nhUi^^ 
ik^ik^lkU=sy>.  ^  ^ 


ik=\  jk-\  lk=\ 


^0  +7’,  V  rjGJ  W,  dy 
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Nh^ 

1  X--  rr^  rj^  P'=%  mhL-^+mhLjj^  +  mhLn^  tihL-,^  +  nhLj^  +  nhL,i^ 

■32.  2-  l.^ikf^k^j]y=sy,^  y 

ik=l  jk=l  lk=l 


■  T-o.-^T-; 


,/)(wf[G.]W',)rfy] 


Let  the  r,/th  term  of  the  matrices  Wq  and  Wj  be 


WAr,i)  =  Wo{r,i)  (Sly  +  S2) 


mW0{r,i)  nWO(r,i) 


Ti/  /  -A  M/  /  -x  /'Cl  ,  co^ '”^1  n»Vl(r,/) 

Wj  (r,  0  =  (r,  0  (51}'  +  -^2)  >' 


(B.178) 


(B.179) 


respectively.  Substituting  the  spar  line  equation  from  Eq.  A.57  for  all  x  in  Eq.  B.177,  the 
mass  matrix  term  in  the  /th  row  and  yth  column  for  the  jlth  spar  web  layer  may  be 
expressed  as 


(r,p) 


r=lp=l 


^yR-syt 


-Nh^ 

2  '■)  ivo  (p.i)  ^1/)  00 + i2)  "''‘/"Vy 

-ik=\ 


NhyNhy 


4t  s«XC^(^<)/,+ 


ik=l  jk=l 


+W,(r,i)WoipJ)  (51y  +  /S2)'"%”^^  jdy 


Nhi^  N\ 


-5S  + (51y  +  52)'"“V 


ik=\  jk=l 


mhl  nL2 


mL3  I^^  (B.180) 

+Wi  (r,  0  Wo  (pj)  iS^y  +  S2)  dy 
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where 


mUl  =  mhUi^.  +  mWO  (r,  i)  +  mWQ  (j>,j) 
nUl  =  nhU.j^  +  nWO  (r,  i)  +  nWO  {p,j) 
mLl  -  mhL.i^  +  mWO  (r,  i)  +  mWO  (p,j) 
nL\  =  nhL.i^  +  nWO  (r,  i)  +  nWO  (p,j) 
mU2  =  mhU^i^  + mhUji^  + mWO  ir,i)  +mWl  (pj) 
nU2  =  nhU.i^  +  nhUji^  +  nWO(r,i)  +nWlip,j) 
mL2  =  mhL.j^  +  mhLjj^  +  mWQ  (r,  i)  +  mW\  (,pj) 
nL2  =  nhL^i^  +  nhLji^  +  nWO  (r,  i)  +nWl  (p,j) 
mU3  =  mhU.f^  +  mhUji^  +  mWlirJ)  +mWOip,j) 
nU3  =  nhU.i^  +  nhUji^  +  nWl  (r,  i)  +nWO{p,j) 
mL3  =  mhL.j^  +  mhLji^  +  mWl  (r,  i)  +  mWO  (p,j) 
nL3  =  nhL.j^  +  +  nWl  (r,  i)  +  nWO  (pJ) 

mU4  =  mhU-j^  + mhUji^  + mhUii^  + mWl  (r,  i)  +mWl  (p,j) 
nU4  =  nhU^j^  +  nhU.,^-\-nhUi^  +  nW\{r,i)  +nW\  {p,j) 


(B.181) 

(B.182) 

(B.183) 

(B.184) 

(B.185) 

(B.186) 

(B.187) 

(B.188) 

(B.189) 

(B.190) 

(B.191) 

(B.192) 

(B.193) 

(B.194) 
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mLA  =  nihL-i^  +  rnhLjj^  +  nihLij^  +  mWl  (r,  i)  +  mW\  ip,j) 


(B.195) 


nLA  -  nhL.i^  +  nhLji^  +  nhLii^  +  hWI  ir,i)  +nW\  {p,j) 


(B.196) 


Noting  that  Eq.  B.180  is  composed  of  linear  combinations  of  the  integral  family  Isp{m,n) 
defined  in  Eq.  A.70,  the  final  expression  for  the  ijth  term  of  [K^blji  is 

5  5 

r=l p=\ 


X  (r,  0  Wo  (PJ)  [  To/sp  +  Wsp  +  1) 


2  (r,  i)  Wo  (pJ)  I  TqJsp  (mLl,  nil)  +  (mLl,  nil  +  1) 


Nh^,  Nhfj 


2^ 

ik=l  jk-l 


S  I  ^ik^jk[  ^0  0  iP^j)  To.Jsp  (mU2,nU2) 


+Wo  (r,  0  Wi  (pJ)  fi  hp  imU2,nU2  +  1) 


NhyNhy 


ik=\  jk=^\ 


+5  E  Z  ”'o  (P.j)  \lsp  ("■ro,  nt/3) 


tVi',  (r,  i)  Wo (pJ)  fi  Kp (rnUi.nm  +  1) 


SEww'o  (r,  i)  Wj  {p,j)  TqI^p  (mL2,  nL2) 


ik=\  jk=\ 


Nh^  Nh^ 


+Wo  (r,  i)  Wi  (pJ)  f,  hp  (mL2,  nL2  +  1) 


(B.197) 


A  X  Z  '■)  ^0  (A 7)  ^o/s/>  nL3) 


ik=\  jk=l 


+Wi  (r,  0  Wo  (pJ)  ^i./cp  (mL3,  nL3  +  1) 
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NhfjNhyNh^j 

ik=l  jk=i  lk=l 

+fiJ^p(mU4,nU4+l)  j 

Nh^  Nh^  Nh^ 

-|l  X  ('■--■) 'Vi  {p,j)[to/,,{mL4.nL4)) 

ik=ljk=:\lk=l 

+f  j  Igp  {mL4,  nL4  +  1 )  ^  ] 

The  mass  matrix  for  all  the  spar  webs  combined  [K^b^tot  is  obtained  by  summing  the  con¬ 
tributions  of  all  the  spar  webs  given  by  Eq.  B.197.  This  matrix  is  then  appropriately 
merged  into  the  global  stiffness  matrix 


B.13  Mass  Matrix  Contributions  of  a  Rib  Web  Layer 

A  rib  web  is  positioned  in  the  vertical  plane  between  parallel  rib  caps  on  the  upper 
and  lower  wing  surfaces.  This  plane  is  the  x-z  plane  located  aty=ys/B-  The  contribution  to 
the  kinetic  energy  7  of  an  infinitesimal  element  dx  by  dz  of  thej/th  rib  web  layer  is 


dTji  =  2pji^ji(x,z) 


(B.198) 


where  is  the  constant  material  density  of  the  web  layer,  the  velocity  vector  is  the  time 
derivative  of  the  displacement  vector  defined  in  Eq.  B.16,  and  is  the  thickness  of  the 
layer.  The  layer  thickness  may  be  expressed  as  a  linear  function  of  x  only  given  by 


=  ^0  +Ti  X 


(B.199) 


Substituting  Eq.  B.16  into  Eq.  B.198  gives 


(B.200) 


Integrating  Eq.  B.200  over  the  area  of  the  web  layer  now  gives 
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+Z|  ■sX 


+  4<s,)] 


(B.201) 


{q}  dxdz 


where  the  limits  of  z  integration  are  the  depth  distributions  of  the  upper  and  lower  rib  caps 
given  respectively  by  Eqs.  B.133  and  B.134.  Equivalent  to  Eq.  B.200  is  the  expression 


+z(iX)  +  z^(5[5i)]  {q}dzdx 


(B.202) 


where  co  is  the  natural  vibrational  frequency  of  the  web  layer.  Application  here  of  the 
variational  extremum  condition  from  Eq.  1.10  yields  the  mass  matrix  for  the  jlth  rib  web 
layer  given  as 


Performing  the  z  integration  analytically  using  Eq.  B.137  gives 


(B.203) 


{X,  y)  -  hi  U.  y) 
+^(^hl,(x,y)  -A’(j:,y) 


Substituting  Eq.  B.199  for  the  layer  thickness  and  Eqs.  B.133  and  B.134  for  the  spar  cap 
depth  distributions  gives 


^^rwh^jl  ~  Pjl 


■Nh,. 


U  cx=rx,  mhU:.  nhU-J  ^ 


y 


dk=l 


Tq  +T-,x 

^jl  b'/ 


dx 


Nh, 


s«y:: 


L  rjc=rjc^  nhL, 


X 


y 


ik~l 


i 


1 


\dx 
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Nhy  Nhy 


2  ^  ^  jkjx=rxp 


U,.Urx=rx^  mhU-^i^  +  mhUji^  nhU^i^  +  nhUj^ 


Tq.^  +  7’j^.xJ[^|^5qSj  j  +  [^5j5q  j  jdx 


ik=l  jk=l 


NK 


V  V  f'*” 

2i  2i  ^ik^jk]x= 


rx^  mhLii^+mhLji^  nhli^  +  nhLjJ  ^  ^ 

y 


5o  I  dx 


ik=l  jk=\ 


NhyNhyNh^ 


5 1  s 


U^.U^Mrx=rx^  mhU.^  +  mhUji^  +  mhUi^  nhU +  nhU +  nhU ^ 


ik=i  jk~l  lk=l 


Nhj^  Nh^  Nh^ 


(B.205) 


lx-'  K'  ZJ^  ZJ^  zr^  F="a 

2  ^  ^  /^Jj:=r;Kr'^  ^ 


fA:=l  y^=l  /^=1 


Jo  +  r,  ;c  S]S,  dx 

'^ii  •';;  All/ 


Since  y  is  equal  to  a  constant  ynjB  for  the  rib  line  defined  by  the  upper  and  lower  rib  caps, 
the  mass  matrix  term  in  the  /th  row  and  yth  column  for  the  yVth  spar  web  layer  may  be 
expressed  as 

3 

TT^  rjc— A  ^  ^  ntUl  nUl  , 

P/(E 

r=l  -ik=l 


1-0,,- vr'o- 


NhyNh^ 


+  5  ^  ^  ^ik^jk]x=rxX  \  ^RIB  +  ^  yRlB 


ik=l  jk=l 


Nhi  Nh^  . 

1  X-'  ZJ^ZX^  F=’'^a{  ^  ^  „Y  mL2  nL2  mL3  nL3 


■2Z  ll^ik^^jkjx=rx,y^^j>  ■  "1/ 

ik=l  jk~l 

NhjjNhyNhfj 


^0,;  +  V  ^RIB  ^  yRIB 


1 1 2 

ik=ljk=llk~l 


(B.206) 


Nh^  Nh^  Nh^ 


ik=l  jk=\  lk=l 
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where  the  x  and  y  powers  are  defined  in  Eqs.  B.141-B.156.  Noting  that  Eq.  B.206  is  com¬ 
posed  of  linear  combinations  of  the  integral  family  /^(m,n)  defined  in  Eq.  A.  102,  the  final 
expression  for  the  ijth  term  of  [M^i,]ji  is 


3  rNhu 


r=l  Lilc=l 


Py/X  X^4  +  l,nUl) 


X^/  +Ti  hg(mLl  +  I,  nil) 


Nh^  Nh^ 


+5  E  E  V*s  *  ^h/nB  <"■“ + "^2) 


2^  ^ 

ik-l  jk-l 


Nhy  Nhy 


X  X  {mm,  nU3)  +  (mU3  +  1,  nU3) 

ik-l  jk=l 


Nhi  Nh^ 

”2  S  ^  ^ RB  (^^2,  nL2)  4-  {mL2  +  1,  nL2)  1 

ik=l  jk=l 
Nhi  Nh^ 

ik=\  jk=l 

Nh^j  Nhy  Nhu 

X  X  X  ^^ik^^k^lijo^RB  nf/4)  +  {mUA  +  1,  nUA)  1 

ik=\  jk=i  lk=\ 

Nh^  N\  Nh^  1 

1  E  E  E  V«»  +  K’kb  *'‘.nIA)) 

ik=ijk=llk=l 


(B.207) 


The  mass  matrix  for  all  the  rib  webs  combined  is  obtained  by  summing  the  contri¬ 

butions  of  all  the  rib  webs  given  by  Eq.  B.207.  This  matrix  is  then  appropriately  merged 
into  the  global  mass  matrix 
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B.14  Stiffness  Matrix  Contributions  of  a  Rib  Web  Layer 

A  rib  web  layer  in  the  ;c-z  plane  is  treated  as  a  plane  stress  panel  where  the  only  strains 
of  importance  are  and  From  the  assumptions  stated  in  Section  B.2,  may  be 

neglected.  Using  Eqs.  B.4,  B.7,  B.9,  B.l  1,  and  B.13  let  us  define 


where 


{  [FQ(x,y)]  +z[Fj  (jc,y)]  }  {q{t)} 


(B.208) 


(B.209) 

(B.210) 

(B.211) 


The  matrices  Fq  and  Fj  are  2  x  {Nq^+Nq-i+Nq^)  in  dimension  and  contain  polynomial  terms 
oftheform  C(r,0^  y  ^  ^  and  C(r,0:t  ’y  respectively.  They 

are  partitioned  into  subvectors  of  dimension  1  x  Nq^  (5=1,3,5). 

Now  the  contribution  to  the  potential  or  strain  energy  U  of  an  infinitesimal  element  dx 
by  dz  of  the  y7th  rib  web  layer  is 


dUji  -  (x,  z) 


dxdz 


(B.212) 


where  [Sly/  is  the y7th  layer’s  constitutive  matrix,  and  tji{x,z)  is  the  thickness  of  the  layer. 
The  constitutive  matrix  is  defined  identically  to  that  of  the  spar  web  in  Eqs.  B.166-B.169 
where  the  T|  axis  is  defined  to  be  the  x  axis.  The  layer  thickness  may  be  expressed  as  a 
function  of  x  only  as  given  in  Eq.  199.  Substitution  of  Eq.  B.208  into  Eq.  B.212  gives 

dUj,  =  kw  {9}’'K[01;,Fo  +  2('^[2]a) 

(  t  \  2(  T  'll  - 

+4^;  [2]^,F„J  +  Z  (f;  [S]„F,J]  {qidxdz 


143 


Integrating  this  over  the  area  of  the  web  layer  now  gives 


+z  f[ [S]^/F J  +  z\ f\  [Q] {~q}  dzdx 


(B.214) 


where  the  limits  of  z  integration  are  the  depth  distributions  of  the  upper  and  lower  rib  caps 
given  in  Eqs.  B.133  and  B.134.  Application  here  of  the  variational  extremum  condition 
from  Eq.  1.10  yields  the  stiffness  matrix  for  the  jlth  rib  web  layer  given  as 


+z[  f[  [  m  jiF^  1  +  z1  F;  [  0]  jiF, )  ]  dzdx 


(B.215) 


Performing  the  z  integration  analytically  using  Eq.  B.137  gives 


//  “  \x=rxhl^  I-  ^  ^0  ^21 


^  f  n  9  It  \  f  t 

hi  ix,  y)  -  hi  (X,  y)  F,  + 1  F,  [^jiF^ 


+^[hlix,y)  -hlix,y)j\^Fl[U]jiF^ j^dx 


(B.216) 


Substituting  Eq.  B.199  for  the  layer  thickness  and  Eqs.  B.133  and  B.134  for  the  rib  cap 
depth  distributions  gives 

\  /  \ 


■KIZ 


L  (x=rx.  mhL;^  nhL^J  .  .  V  r-?’  i 


y  “  To  +F,  X  f;[S]^zFJ^ 


Nh^Nhij 


1  V  Tj^rr^F=^^A  mhU-^  +  mhUj!^  nhU.^  +  nhU^^ 

’2  2^  ^ik^jk]x=rxp^  ^ 

ik=l  jk=\ 


f.  +  Kx  fJ  [  m  jiF,  +  f[  [  S]  jiF, )  )dx 
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NK  NK 


rx^  mhLii^  +  mhLji^  nhL-^  +  nhLj^ 


ik=\  jk-\ 


+  vjll^o  [21,7^ J  +  [f] m]jiF, ] )dx 


(B.217) 


NhyNh^Nhy 


j^U  j^U  j^U  rx=rxj^  mhU +  mhUjf^  +  mhU nhU-f^  + nhUjf^  + nhUn^ 
^  ^  ik  jk  lkjx=rXp^  ^ 


ik=l  jk=llk=\ 


Nh^  Nh^  Nh^ 


1  1  7  T  rx=rx^  +  mhLj}^  +  mhLi^  nhL-f^  +  nhLj,^  +  nhLi^ 

3  ^  ^  ^  ik  jk  lk]x=rxp^  ^ 

ik=l  jk-\  Ik-l 


•rO,+  VA^lt2]7/^l 


Let  the  r,ith  term  of  the  matrices  Fq  and  F,  be 


„  ,  ~  mFO{r,i)  nF0(r,i) 

FQ{r,i)  =  FQ{r,i)x  > 


Fj(r,0 


(B.218) 


(B.219) 


respectively.  Since  y  is  equal  to  a  constant  along  the  spar  line  defined  by  the  upper  and 
lower  rib  caps,  the  mass  matrix  term  in  the  ith  row  and  yth  column  for  the  y/th  spar  web 
layer  may  be  expressed  as 


r=l  p=l 


yNK  ,  X 

X  Hlh  (n  0  h  (pj)  IZZiK  * 

^ik=l 


F  /  X 

■  X  H^Foir,  i)  h  (pj) 

ik=l 
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Nh^  Nhy 


4.1  V  V 


rx=r 

Jx=r 


ik=l  jk-l 


Nh^  Nh^  . 

-2  I  Z  F,  (pJ)x"“'M% 


^  ^  ~  ,  mill  nUl 

Tq.^  +  T 1 .  X J Fq  ( r,  i)  Fi  (p,  j)  X 

~  .  -  ,  mm  nm  ^  . 

+F^{r,i)FQ{p,j)x  y^ig  jdx 


ik=l  jk=\ 


•^  c  /'  •^  (B.220) 

+Fiir,i)FQip,j)x  yj^jg  jdx 


+ 


Nh^Nh^Nh^ 

5  z  z  z  (-•.'■)  '"i  (P’i)  izzi  fo, + 

ik=ljk=llk=l 

Nhi^  Nhi^  N\ 

I  z  z  z  (-■.  -■)  h  (p.i)  jz::i\ *  v>”“ 

ik=ljk=llk=l 


where 


mill  =  mhU-i^  +  mFO  (r,  i)  +  mFO  (p,j) 

(B.221) 

nUl  -  nhU-j^  +  nFO  (r,i)  +nFO(p,j) 

(B.222) 

mLl  =  mhL-j^  +  mFO  (r,  i)  +  mFO  ip,j) 

(B.223) 

nLl  =  nhL-j^  +  nFO  {r,  i)  +  nFO  (p,j) 

(B.224) 

mU2  =  mhU./^  +  mhUji^  +  mFO  (r,  /)  +  mFl  (p,j) 

(B.225) 

nU2  =  nhU.i^  +  nhUji^  +  nFO  (r,  i)  +nFl(p,j) 

(B.226) 

mL2  =  mhL.i^  + mhLji^  + mFO  (r,  i)  +  mFl  ip,j) 

(B.221) 

nL2  =  nhL-i^  +  nhLjj^  +  nFO  (r,  i)  +  nFl  {,p,j) 

(B.228) 

mU3  =  fnhU-i^  +  mhUji^  +  mFl  (r,  i)  +  mFO  (p,j) 

(B.229) 

nU3  =  nhU.i^  +  nhUji^  +  nFl  (r,  i)  +nFO(p,j) 

(B.230) 
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mL3  =  mhL,.+mhL..+mF[ir,i)  +mFOip,j) 

Ik  jk  (B.231) 

nL3  =  nhL..+nhL.,+nFl(r,i)  +nFO(p,j) 

(B.232) 

m(/4  =  mhU;.  +  mhU-.  +  mhU,.  +  mFl  (r,  i)  +  mF\  (pj) 

Ik  Jk  Ik  (B.233) 

nU4  -  nhU..+nhU.,+nhU,.+nF\{r,i)  +nF\{p,j) 

Ik  Jk  Ik  (B.234) 

mL4  =  mhL..  +  mhL.,  +  inhL,.  +  mFl  (r,  i)  +  mFl  (pJ) 

Ik  Jk  Ik  (B.235) 

nL4  =  nhL..  +  nhL,,  +  nhL,.  +  nFl  (r,  i)  +  nFl  (pJ) 

Ik  Jk  Ik  (B.236) 

Noting  that  Eq.  B.220  is  composed  of  linear  combinations  of  the  integral  family  Ii^g{m,n) 

defined  in  Eq.  A.  102,  the  final  expression  for  the  ijth  term  of  [K^y^blji  is 


ft-jj 

Z  E  ('’>-'■)  \'kb  nU2) 

ik=l  jk=l 

+Fo  (r,  0  Fj  (pJ)  {mU2  +  1,  nt/2)  j 

Nh^  Nhy 

+5  z  z  Vs" 

jk=l 

+Fi  (r,  0  Fq  (pJ)  TiJjfg  imU3  +  1,  nU3)  ^ 
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Nh^  Nh^ 


Fq  (r,  i)  F 1  {p,  j)  Tq {mil,  nil) 

\ 

ik=l  jk~l 

+FQ{r,i)F^{p,j)f^J^g{mLl  +  \,nL2)  J 

Nh,  Nh^ 

(B.237) 

Fj  (r,  /)  Fq  (pJ)  Fq imL3,  nL3) 

ik^l  jk=l 

+Fi  (r,  i)  Fq  ipj)  (mL3  +  1,  «L3)  1 

NhjjNhfjNh^j 

*{111  "Ifilifllifi  (>■.  0  ^1  (P’J>  [  V«s 

ik-\ jk-\  lk=\ 

+ri./^g(mf/4+  l,nf/4)  j 

Nh^  Nh^  Nh^ 

X  \ 

iIXX4 

('•’  0  h  (pJ)[\Jrb  ("^^4,  nL4)  J 

ik=l  jk=llk=l 

h-Ti  /^g  (mL4  +  1,  nL4)  j  J 

The  stiffness  matrix  for  all  the  rib  webs  combined  is  obtained  by  summing  the 

contributions  of  all  the  layers  of  all  the  rib  webs  given  by  Eq.  B.237.  This  matrix  is  then 

appropriately  merged  into  the  global  stiffness  matrix 


Appendix  C: 

Fundamental  Integrals  and  Tables 


C.1  Introduction 

Integrating  energy  contributions  over  the  geometry  of  wing  components  has  led  to  the 
definition  of  three  integral  families:  IxR{tn,n),  Ispini,n),  and  lRs{fn,n).  This  appendix  dis¬ 
cusses  the  analytical  solutions  to  these  integrals  based  on  the  wing  model’s  shape  design 
variables.  The  discussion  starts  with  the  development  of  a  family  of  solutions  to  integrals 
over  the  area  of  a  trapezoidal  panel.  Next  is  the  development  of  a  family  of  solutions  to 
integrals  over  the  length  of  a  spar.  Last  is  the  development  of  a  family  of  solutions  to  inte¬ 
grals  over  the  length  of  a  rib.  From  these  solutions,  solution  tables  of  each  integral  family 
may  be  constructed.  These  tables  need  only  be  calculated  once  for  a  particular  wing 
geometry  and  may  then  be  referenced  multiple  times  in  mass  and  stiffness  matrix  calcula¬ 
tions. 


C.2  Area  Integrals 

Derivation  of  mass  and  stiffness  matrix  equations  for  trapezoidal  skin  panels  has 
resulted  in  the  definition  of  a  family  of  integrals  given  by 

=  llx^y^dxdy 

where  the  trapezoidal  area  of  x-y  integration  is  defined  by  the  vertices  {xpi,yi)y  {XpRyy^, 
(xAR^yRl  U/U-Jr)  as  shown  in  Figure  1.3.  The  left  and  right  sides  of  the  trapezoid  are  paral¬ 
lel  to  the  jc  axis  and  located  at  y  coordinates  and  yg  respectively.  Equations  for  the  front 
and  aft  lines  of  the  trapezoid  may  be  written  as 


Xj,iy)  =  Fly  +  F2 


(C.2) 


(y)  =  A\y  +  A2 


(C.3) 
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respectively,  where 


^FR  ^FL 

yR-yi 


(C.4) 


P2  =  ^FRyi 

yR-yi  (C.5) 


A1  = 


^AR  ^AL 

yR-yi 


A2  = 


^AiyR-^ARyL 

yR-yL 


The  integral  may  now  be  written  as 

V=y«  nK=Xp{y)  m 


(C.6) 


(C.7) 


(C.8) 


Using  Eqs.  C.2  and  C.3  gives  the  relationship 

=  ^_^[(Aly+A2)'”^^-  (Fly  +  F2)'”^T 

ix=Xp(y)  m+1 

Substituting  Eq.  C.9  into  Eq.  C.8  gives 

Vs*'".'’)  =  + 

+j^:yy(^iy+’^2)"‘^Vy] 

Let  us  define  two  secondary  integrals  7^  and  Ip  where 

V  ('■.*)  =  fZl‘/(Aly+A2)‘dy 


and 


(C.9) 


(C.lO) 


(C.ll) 


(C.12) 
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Recursion  formulas  may  be  used  to  obtain  a  closed  form  evaluation  of  these  secondary 
integrals  (SP68,LIV90).  They  are  evaluated  using  the  formulas 

J/(Aly +A2)''^f>'= 

(C.13) 

y'*'  +  +  ^^1/ (Aly+A2)-'dy 

r  +  5+1  r  +  s'+W 


and 


j/iFly  +  Fiydy-- 


(C.14) 


/^^Fly  +  F2)^  ^  sF2 
r  +  5:+  1  r  +  s  + 


-j/iFly  +  Fiy-'dy 


Substitution  of  Eqs.  C.ll  and  C.12  into  Eq.  C.IO  yields  a  final  expression 

Irj,„{m,n)  =  — r  [/.  (n,  m+ 1) m  +  1) ] 

TR  m+l  ^  ^  (C.15) 

Tables  for  and  Ip  may  be  constructed  for  ascending  x  and  y  powers  using  the  formulas  in 
Eqs.  C.13  and  C.14.  From  these  a  final  table  may  then  be  constructed  for  the  trapezoidal 
area  integrals  using  Eq.  C.15  (LIV90). 


C.3  Spar  Line  Integrals 


Derivation  of  mass  and  stiffness  matrix  equations  for  spars  and  spar  webs  has  resulted 
in  the  definition  of  a  family  of  integrals  given  by 


(m,  n)  =  - r  (51y  +  SI)  ^y^dy 

syj^  -  sy^Jy-^yi 


(C.16) 


where  the  spar  line  endpoints  ai'e  defined  by  the  coordinates  and  isxp,syii)  as 

shown  in  Figure  1.4,  and  L  is  the  length  of  the  spar.  This  integral  results  form  substituting 
for  X  the  spar  line  equation  given  by 
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jc(y)  =  51}’  +  52 


(C.17) 


where 


51  = 


SXj^  -  SXj_^ 

syR-syi 


(C.18) 


52  = 


sXiSyR-sxRsyi 

^yR-^yi 


(C.19) 


Define 


4 


(C.20) 


which  may  be  evaluated  in  closed  form  using  the  recursion  formula 


r  +  s+  1  r  +  ^  +  IJ 


(C.21) 


Substitution  of  Eq.  C.20  into  Eq.  C.16  yields  a  final  expression 


lRp{m,n)  - 


L 


^yR-^yi 


(m,  n) 


(C.22) 


A  table  of  spar  line  integrals  for  ascending  x  and  y  powers  may  be  prepared  using  the  for¬ 
mula  from  Eq.  C.21  combined  with  Eq.  C.22. 


C.4  Rib  Line  Integrals 


Derivation  of  mass  and  stiffness  matrix  equations  for  ribs  and  rib  webs  has  resulted  in 
the  definition  of  a  family  of  integrals  given  by 


/^g(m,  n)  -  yR[Blx=rXp^’^^^ 


(C.23) 


where  the  rib  line  endpoints  are  defined  by  the  coordinates  (rxFynjB)  and  (rx^ym)  as 
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shown  in  Figure  1.4.  The  coordinates  rxp  and  are  determined  from  the  rib  shape  design 

variables  using 


rXp  —  Fiy  +  F2 


rx^  —  Aiyj^jg+A2 


(C.24) 

(C.25) 


where 


FI  = 


rXpR-rXpL 


(C.26) 


F2  = 


rXpjyR-rXpRryL 

^yR-^yi 


(C.27) 


A1 


^^AR-^^AL 

ryR-ryi 


(C.28) 


A2  = 


r^ALn'R-fXARryi 

ryR-ryp 


(C.29) 


The  rib  integral  is  easily  evaluated  and  given  by 


.  .  .  yRlB  r  m+l  m+l-l 


(C.30) 


This  equation  may  be  used  to  construct  a  table  of  rib  integrals  for  ascending  x  and  y  pow¬ 


ers. 


Appendix  D: 
Test  Case  Data 


D.1  Overview 

Contained  in  this  appendix  is  the  structural  data  for  the  wing  test  cases.  Data  for  the 
Turner-Martin- Weikel  wing  is  presented  first.  Given  next  is  the  data  for  the  HSCT  wing. 


D.2  Turner-Martin-Welkel  Wing 

Figure  D.l  shows  a  finite  element  model  of  the  Tumer-Martin-Weikel  (T-M-W)  wing 
(LIV194).  Figure  D.2  illustrates  details  of  the  construction  of  the  wing  (TMW64).  Depth, 
skin  thickness,  spar  web  thickness,  and  spar  cap  area  are  all  constant.  The  wing  is  made 
totally  from  aluminum  having  the  following  material  properties: 

E  =  10.5E6  ;  v  =  0.3  ;  p  =  0.1  lb/in^ 


Figure  D.  1:  FEM  Model  of  the  T-M-W  Wing 


TYPICAl  CHOftOWISE  SCaiON  THROUGH  MOOCL 


Figure  D.2:  T-M-W  Wing  Construction 
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D.3  HSCTWIng 

Figure  D.3  shows  the  ELFINI  finite  element  model  of  the  Boeing  HSCT  candidate 
wing.  The  wing  is  symmetric  about  the  mid-plane  with  the  depth  distribution  for  the 
upper  skin  given  in  inches  by 

h  =  25.0  -  0.03097243y 

The  wing  skins  are  composed  of  composite  material  layers  with  the  following  material 
properties; 

=  36.2E6  psi;  £22  =  1-4E6  psi;  v  =  0.29;  =  .666E6  psi;  p  =  0.1  Ib/in^ 

The  inboard  upper  and  lower  wing  skins  are  each  made  of  20  sets  of  (0/90/45/-45)  lami¬ 
nate  where  ply  orientation  is  referenced  to  the  inboard  spar  lines.  The  same  construction 
is  used  for  the  outer  wing  where  the  plies  are  oriented  with  respect  to  the  outer  wing  spar 
lines.  The  thickness  of  each  ply  is  0.0037  in.  Spar  and  rib  webs  are  made  of  the  same 
material  using  4  sets  of  (0/90/45/-45)  laminate  where  orientation  is  referenced  from  the 
mid-plane  (LSB93). 


Figure  D.3:  ELFINI  Model  of  HSCT  Wing 
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Spar  and  rib  caps  have  a  constant  area  of  0.267  in^.  They  are  made  of  composite 
material  with  the  following  properties: 

E  =  I4.0E6  psi;  v  =  0.46;  p  =  0.1  Ib/in^ 

Table  D.  1  lists  the  vertical  spring  coefficients  used  at  points  along  the  root  to  represent 
fuselage  stiffness  along  the  side  of  body  rib  (LSB93).  Linear  springs  with  a  coefficient  of 
2.0E10  Ib/in  are  used  along  the  root  to  impose  zero  Vo  displacement.  Rotational  springs 
with  coefficients  of  1.0E16  in-lb/rad  along  the  root  aty=83  in.  and  1.25E16  in-lb/rad  along 
the  carry  thi'ough  root  are  used  to  impose  zero  tj/j,  displacement.  One  spring  with  a  coeffi¬ 
cient  of  I.OEIO  is  used  at  the  rear  inboard  point  of  the  wing  carry  through  to  resist  Mq  dis¬ 
placement 


Table  D.l:  Root  Springs  Representing  Fuselage 
Stiffness  Along  Side  of  Body  Rib 


Spring  # 

X  (in) 

%  (ib/in) 

Spring  # 

X  (in) 

%  (Ib/in) 

1 

1393.27 

3257 

17 

1939.82 

20127 

2 

1428.43 

3572 

18 

1975.02 

23397 

3 

1463.58 

3917 

19 

2010.22 

30295 

4 

1498.74 

4351 

20 

2045.22 

39460 

5 

1533.90 

4850 

21 

2080.22 

48539 

6 

1569.06 

5396 

22 

2115.22 

59351 

7 

1604.22 

6035 

23 

2150.22 

73546 

8 

1631.22 

6589 

24 

2177.95 

85157 

9 

1658.22 

7202 

25 

2205.68 

108019 

10 

1693.42 

8137 

26 

2238.72 

120000 

11 

1728.62 

9525 

27 

2269.05 

171571 

12 

1763.82 

10359 

28 

2299.38 

173602 

13 

1799.02 

11650 

29 

2329.71 

169070 

14 

1834.22 

13267 

30 

2360.04 

162117 

15 

1869.42 

15174 

31 

2390.37 

149298 

16 

1904.62 

17423 

32 

2420.70 

100000 

Appendix  E: 

CONNECT  Program  Information 


E.1  Overview 

This  appendix  briefly  describes  the  code  developed  and  used  in  this  work.  Program 
subroutines  are  listed  with  a  description  of  their  purpose.  The  program  structure  is  out¬ 
lined  with  a  description  of  the  main  program  variables. 


E.2  Program  Subroutines 

CONNECT,  the  main  program,  is  written  in  Fortran.  Following  is  a  list  of  the  sup¬ 
porting  subroutines  grouped  according  to  their  purpose. 


Matrix  Preparation: 
MTXSKN  - 
MTXRIB 
MTXSWB  - 
MTXRWB  - 
MTXMAS  - 
MTXCORE  - 


prepares  matrices  for  FSDPT  skin  and  spar  stiffness  calculations 
prepares  matrices  for  FSDPT  rib  stiffness  calculations 
prepares  matrices  for  FSDPT  spar  web  calculations 
prepares  matrices  for  FSDPT  rib  web  calculations 
prepares  displacement  matrices  for  FSDPT  mass  calculations 
prepares  matrices  for  FSDPT  equivalent  core  calculations 


Stiffness  Matrix  Evaluation: 

SKIN  -  evalutes  stiffness  matrix  contributions  for  a  single  FSDPT  skin  layer 

SPRCAP  -  evaluates  stiffness  matrix  contributions  for  a  single  FSDPT  spar 

RIBCAP  -  evaluates  stiffness  matrix  contributions  for  a  single  FSDPT  rib 

SPRWEB  -  evaluates  stiffness  matrix  contributions  for  a  single  FSDPT  spar  web 

layer 

RIBWEB  -  evaluates  stiffness  matrix  contributions  for  a  single  FSDPT  spar  web 


layer 

SKINO  -  evaluates  stiffness  matrix  contributions  for  a  single  CPT  skin  layer 

CORE  -  evaluates  stiffness  matrix  contiibutions  for  an  equivalent  core 


SPNGIG  -  evaluates  stiffness  matrix  contributions  for  a  FSDPT  zone  attached  to 
ground  at  a  point  via  springs 

SPNGl  1  -  evaluates  stiffness  matrix  contributions  for  the  spring  connection 

of  2  FSDPT  zones  at  a  point 

SPNGIO  -  evaluates  the  stiffness  matrix  contributions  for  the  spring  connection 
of  a  CPT  zone  to  a  FSDPT  zone  at  a  point 

Mass  Matrix  Evaluation: 

MSKIN  -  evaluates  mass  matrix  contributions  for  a  single  FSDPT  skin  layer 

MSPRCAP  -  evaluates  mass  matrix  contiibutions  for  a  single  FSDPT  spar 

MRIBCAP  -  evaluates  mass  matrix  contributions  for  a  single  FSDPT  rib  cap 

MWEB  -  evaluates  mass  matiix  contributions  for  a  single  FSDPT  spar  or  rib 

\veb  layer 

MCONC  -  evaluates  mass  matrix  contribution  for  a  single  FSDPT  concentrated 
mass 

MSKINO  -  evaluates  mass  matrix  contributions  for  a  single  CPT  skin  layer 
Integral  Table  Generation: 

TABSRF  -  generates  integral  table  for  a  single  panel  geometry 
GMTRY  -  used  in  TABSRF  to  calculate  panel  geometry  parameters 
INTGSP  -  generates  integral  table  for  a  single  spar 

INTGRB  -  generates  integral  table  for  a  single  rib 

Skyline  Linear  Solver: 

SKYFAC  -  factorizes  global  stiffness  matrix 
SKYFAC  -  solves  for  static  generalized  displacements 


Eigensolution: 

EIGZF 


-  solves  eigenvalue  problem  for  natural  frequencies  and  mode  shapes 


E.3  Program  Structure 

Phase  1:  Initialization 
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All  input  parameters  are  set  within  the  main  code  and  stored  in  arrays.  No  data  is  read 
from  file.  To  avoid  numerical  ill-conditioning,  all  geometric  parameters  may  be  shifted 
and  scaled  to  fit  into  a  unit  square.  All  other  parameters  with  length  dimensions  must  also 
be  scaled  for  consistency.  Following  is  a  list  and  description  of  array  variables  with  their 
index  descriptors. 


idzone(zone#) 

nql(zone#)...nq5(zone#) 

mq  lx(term#,zone#)... 

. . .  mq5x  (term#, zone#) 
mq  1  y  (term#,  zone#) . . . 

...mq5y(term#,zone#) 

scale(zone#) 

xshft(zone#) 

yshft(zone#) 


-  identifies  CPT  (0)  or  FSDPT  (1)  zone 

-  #  of  displacement  terms  for  Uq,  Vq,  t]/,,  vvq  respectively; 

(  w  from  CPT  zone  is  treated  as  Wq  from  FSDPT  zone) 

-  a:  powers  for  displacement  polynomials 

-  y  powers  for  displacement  polynomials 

-  zone  scaling  factor 

-  zone  X  coordinate  shift  value 

-  zone  y  coordinate  shift  value 


npanel(zone#) 

xyplan(coord#,panel#,zone#) 

pnlh  (p  anel#,zone#) 

nlayr(panel#,zone#) 

ntterm(layer#,panel#,zone#) 

t(term#,layer#,panel#,zone#) 

mtx(trm#,layr#,pnl#,zone#) 

mty(trm#,layr#,pnl#,zone#) 

skbeta(layr#,pnl#,zone#) 


#  of  panels  in  zone 

panel  geometry  coordinates 
panel  depth  distribution  (h  series)  tag 

#  of  skin  layer  in  panel 

#  of  thickness  distribution  terras  in  skin  layer 
layer  thickness  distribution  coefficient 

X  powers  of  layer  thickness  distribution 
y  powers  of  layer  thickness  distribution 
skin  layer  fiber  orientation  from  x  axis 


nhterm(h  series#,zone#) 
h(terra#,h#,zone#) 
mhx(term#,h#,zone#) 
mhy(term#,h#,zone#) 


-  #  of  terms  in  depth  distribution  (h  series) 

-  depth  distribution  (h  series)  coefficient 

-  X  powers  of  depth  distribution 

-  y  powers  of  depth  distribution 


160 


spyl,spyr,spxl, 

spxr(spap#,zone#)  -  spar  geometry  coordinates 

sprh(spar#,zone#)  -  spar  depth  distribution  tag 

sparea(term#,spar#,zone#)  -  spar  cap  area 

rbyl,rbyr,yrib,rbxfl,rbxfr 

rbxal,rbxar(rib#,zone#)  -  rib  geometry  coordinates 

ribh(rib#,zone#)  -  rib  depth  distribution  tag 

rbarea(term#,spar#,zone#)  -  rib  cap  area 

npnts(zone#)  -  #  of  attach  points  in  zone 

pntxyz(xyz#,point#,zone#)  -  attach  point  geometry 

pnth(point#,zone#)  -  attach  point  depth  distribution  tag 

wbspr(term#,spr  web#, zone#)-  idntifies  spar  caps  associated  with  spar  web 
twbspr(trm#,lyr#,spwb#,zn#)  -  spar  web  layer  tickness  coefficient 
bwbspr(lyr#,spwb#,zn#)  -  fiber  orientation  of  spar  web  layer 

wbrib(term#,rib  web#, zone#)  -  idntifies  rib  caps  associated  with  rib  web 
twbrib(trm#,lyr#,rbwb#,zn#)  -  rib  web  layer  thickness  coefficient 
bwbrib(lyr#,rbwb#,zn#)  -  fiber  orientation  of  rib  web  layer 

cnctn(zone#,zone#)  -  zone  connection  tag:  0-not  connected,  1 -connected 

(ground  connection  when  zone  number  are  equal) 
nspring(zone#,zone#)  -  #  of  spring  used  for  zone  connection 

xattl,xattr,yattl, 

yattr(lor2,zone#,zone#)  -  attach  line  geometry  of  zones  being  connected 
ksprng(trm#,sprng#,zn#,zn#)  -  spring  stiffness  coefficients  corresponding  to 

ncmass(zone#)  -  #  of  concentrated  masses  in  zone 

cmpoint(cmass#,zone#)  -  concentated  mass  point  tag 
cmass(xyz#,cmass#,zone#)  -  concentrated  mass  coefficient 

nlconc(zone#,load  case#)  -  #  of  concentrated  force  loads 
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qconc(ld#,xyz#,zn#,ldcase#)  -  concentrated  force  load  coefficient 
idconc(ld#,xyz#,zn#,ldcase#)  -  concentrated  load  point  tag 
nlpol(panel#,zone#,ldcase#)  -  #  of  distributed  pressure  loads 
qpol(trm#,xyz#,ld#,pnl#, 

zn#,ldcase#)  -  distributed  load  coefficient 
mqpolx(term#)  -  a:  powers  of  distributed  loads 


mqpoly(temi#)  -  y  powers  of  distributes  loads 


noutx(panel#,zone#) 

nouty(panel#,zone#) 

outx(point#,panel#,zone#) 

outy(point#, panel#, zone#) 

yleft(p  anel#,zone#) 

yright(panel#,zone#) 

pnplmt(panel#,zone#) 


-  number  of  panel  output  points  in  x  direction 

-  number  of  panel  output  points  in  y  direction 

-  X  direction  output  grid  coordinates 

-  y  direction  output  grid  coordinates 

-  left  y  coordinate  of  output  grid  on  wing  panel 

-  right  y  coordinate  of  output  grid  on  wing  panel 

-  panel  output  tag:  0-no  print,  1 -print 


noutspr(zone#)  -  number  of  spar  output  points 

outspr(point#,zone#)  -  spar  output  grid  coordinates 

sprprnt(spar#,zone#)  -  spar  output  tag:  0-no  print,  1 -print 


noutrib(zone#)  -  number  of  rib  output  points 

outrib(point#,zone#)  -  rib  output  grid  coordinates 

ribpmt(rib#,zone#)  -  rib  output  tag:  0-no  print,  1 -print 


In  addition  to  these,  material  property  variables  are  also  assigned: 


qskl  1,  qskl2, 
qsk22,  qsk66 
qwbl  1,  qwbl2, 
qwb22,  qw66 
q0skll,q0skl2, 
q0sk22,  q0sk66 
modcap 
modOcap 
densk 


-  FSDPT  skin  constitutive  matrix  values 

-  FSDPT  spar  and  rib  web  constitutive  matrix  values 

-  CPT  skin  constitutive  matrix  values 

-  FSDPT  spar  and  rib  cap  modulus 

-  CPT  spar  and  rib  cap  modulus 

-  FSDPT  skin  density 
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denwb  -  FSDPT  spar  and  rib  web  density 

denOsk  -  CPT  skin  density 

dencap  -  FSDPT  spar  and  rib  cap  density 

denOcap  -  CPT  spar  and  rib  cap  density 

Phase  2:  Calculations 

A)  Integral  Table  Assembly:  The  following  integral  table  arrays  ai'e  assembled. 

aintmn(x  power, y  power, 

panel#,zone#)  -  panel  integral  table 

spintmn(x  power,y  power, 

spai#,zone#)  -  spar  integral  table 

rbintmn(x  power,y  power, 

rib  #,zone#)  -  rib  integral  table 

B)  Stiffness  Matrix  Assembly:  Here  first  subroutines  are  called  to  prepare  the  needed 
strain  matrices.  Each  strain  matrix  from  the  FSDPT  formulation  is  stored  in  3  arrays:  one 
for  coefficients,  one  for  x  powers,  and  one  for  y  powers.  The  arrays  are  defined  as  follows: 

ea, eax,eay(row#,col#,zone#)  -  1st  FSDPT  skin  and  spar  strain  matrix  [7?o] 

eb, ebx,eby(row#,col#,zone#)  -  2nd  FSDPT  skin  and  spar  sti-ain  matrix  [7?i] 

el,elx,ely(row#,col#,zone#)  -  1st  FSDPT  rib  strain  matrix  [To] 
e2,e2x,e2y(row#,col#,zone#)  -  2nd  FSDPT  rib  strain  matrix  [7,] 

esa, esax,esay(rw#,col#,zn#)  -  1st  FSDPT  spar  web  strain  matrix  [Wq] 

esb, esbx,esby(rw#,col#,zn#)  -  2nd  FSDPT  spar  web  strain  matrix  [ITi] 

era, erax,eray(rw#,col#,zn#)  -  1st  FSDPT  rib  web  strain  matrix  [F,,] 

erb, erbx,erby(rw#,col#,zn#)  -  2nd  FSDPT  rib  web  strin  matrix  [FJ 

eec,eecx,eecy(rw#,col#,zn#)  -  Equivalent  core  strain  matrix  [WJ 
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Next  the  subroutines  which  calculate  the  stiffness  contributions  of  each  wing  component 
are  called.  Since  the  stiffness  matrix  is  symmetric,  only  the  upper  triangular  portion  need 
be  stored.  Formatting  of  the  SKYLINE  solver  requires  that  this  upper  triangular  matrix  be 
stored  in  a  vector.  This  requires  a  vector  of  length  Nqgigi,*(Nqgigh+\)l2.  Stiffness  matrix 
terms  are  stored  column  by  column  in  the  vector  array:  stif(term#).  Only  skin  contribu¬ 
tions  are  included  for  a  CPT  zone  in  this  program. 

C)  Loads  Assembly:  Here  the  load  vector  is  assembled  and  stored  in  the  array: 

qload(term#,load  case#) 

D)  Static  Solution:  Here  the  subroutines  for  factorizing  the  stiffness  matrix  and  solving 
the  linear  system  of  equations  are  called.  The  solution  generalized  displacements  are 
stored  in  the  array:  qdisp(term#,load  case#). 

E)  Mass  Matrix  Assembly:  Here  first  a  subroutine  is  called  to  prepare  the  needed  dis¬ 
placement  matrices.  Each  displacement  matrix  from  the  FSDPT  formulation  is  stored  in  3 
arrays  as  the  strain  matrices  are  in  stiffness  matrix  assembly.  These  arrays  are  defined  as 
follows: 

.s0,s0x,s0y(row#, col#, zone#)  -  1st  displacement  matrix  [^q] 
sl,slx,sly(row#,col#,zone#)  -  2nd  displacement  matrix  [SJ 

Next  the  subroutines  which  calculate  the  mass  contributions  of  each  wing  component  are 
called.  The  mass  matrix  is  also  symmetric  and  thus  stored  in  the  same  format  as  the  stiff¬ 
ness  matrix  in  the  vector  array:  mass(term#).  Only  skin  contributions  are  included  for  a 
CPT  zone  in  this  program. 

F)  Eigensolution:  Here  the  eigenvalue  problem  involving  mass  and  .stiffness  is  solved. 
Natural  frequencies  are  sorted  by  increasing  size  and  stored  in  the  array:  eigsort(freq#). 
Corresponding  generalized  modal  displacements  are  stored  in  the  array: 
qmode(term#,mode#). 
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Phase  3:  Output 


Displacement,  stress,  and  mode  shape  output  is  calculated  at  specified  grid  points  and 
stored  in  the  following  arrays: 


dsppnl(disp#,point#,panel#, 

zone#,  load  case#)  - 
dspspr(disp#,point#,spar# 

zone#,load  case#)  - 
dsprib(disp#,point#,panel#, 

zone#,load  case#)  - 
strpnl(stress#,point#,layer#, 
panel#,zone#,load  case#)  - 
strspr(point#,spai#,zone#, 
load  case#) 

strrib(point#,rib#,zone#, 
load  case#) 

mddsppnl(disp#,point#, 

panel#,zone#,mode#) 

mddspspr(disp#,point#, 

spar#,zone#,raode#) 

mddsprib(disp#,point# 

rib#.zone#,mode#) 


panel  displacements 
spar  displacements 
rib  displacements 

panel  stresses  (3  in  x-y  system,  3  in  fiber  axis  system) 
spar  stress 
rib  stress 

panel  mode  shape  displacements 
spar  mode  shape  displacements 
rib  mode  shape  displacements 


Output  is  printed  for  each  panel,  spai’,  and  rib  selected. 


